NUMERICAL REFRACTION AND 
DIFFRACTION FOR COASTAL AREAS 



by 

ROBERT DALTON §MART 

SB, Massachusetts Institute of Technology 
(1957) 



Submitted in partial fulfillment 
of the requirements for the degrees of 
Master of Science in Civil Engineering and Civil Engineer 



at the 

Massachusetts Institute of Technology 
(1967) 



Signature of Author 



Certified by 



Department of Civil Engineering, (May 19, 1967) 



Thesis Supervisor 



Accepted by e « 



Chairman, Departmental Committee on Graduate Students 



fO' _ . 



c , ; / 



a .* 






Ihs 



LIBRARY 

NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CALIF. 93940 

ABSTRACT 



NUMERICAL REFRACTION AND 
DIFFRACTION FOR COASTAL AREAS 



by 

ROBERT DALTON SMART 



Submitted to the Department of Civil Engineering on May 19, 1967, in 
partial fulfillment of the requirements for the degrees of Master of 
Science in Civil Engineering and Civil Engineer. 



The purpose of this study was to develop a computer system which could 
be used by engineers in analysis of water wave behavior in coastal areas. 

A wave ray tracing program was developed which uses three-dimensional 
continuous parabolic interpolation for determination of bottom depths and 
derivatives. Tested on analytical plane beaches this method gave results 
which agreed within 1 per cent of theoretical results. Tested on a 
natural beach, the program agreed favorably with graphical methods and 
numerical methods developed by others* 

For study of diffraction patterns both at ends of semi-infinite 
breakwaters and at breakwater gaps, a program was developed, based on 
Penney-Price methods, which calculates the coefficient of diffraction at 
any point in a diffraction zone. Using formulae derived from the Penney- 
Price methods, a program was developed which calculates directly the di- 
rection of travel of a wave at any point in the diffraction zone, The 
results of this program compare closely with graphical plots of diffraction 
zones which were completed by manual methods. 

These major programs, and a smaller program which locates and traces 
shorelines across an array of water depths, were joined into one system 
through a Problem-Oriented Language. This system reads and works with the 
same raw data available to an engineer in the field performing analysis. 

It will perform tasks and make calculations as directed by the operator, 
then output all information needed by the engineer in his analysis. 

Either analytical or natural beaches, in feet or fathoms, or combina- 
tions thereof may be input to the program. To simulate tidal fluctuations, 
water level changes may be made at any time during analysis. Breakwaters 
may be superimposed on the array of bottom depths; and, if the breakwater 
or its diffraction zones are intersection by a wave during refraction analy- 
sis, all ray values at the point of intersection are calculated. 

Program output includes, shoreline traces, diffraction coefficients 
with direction of travel of diffracted waves, and wave ray traces. Values 
are calculated at given points along the ray including, among other data, 
direction of wave travel, curvature, coefficient of refraction, and the 
shoaling factor. 
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TERMS AND SYMBOLS 



Wherever possible, terms and symbols used in derivatives and text are 
in accordance with the standards of the Council on Wave Research as pre- 
sented by Wiegel [1]\ Hardware and program language requirements prevented 
use of all symbols in actual programs however and the need for additional 
variables may have resulted in notation which is in conflict with the 
standards# Comments within program summaries explain correlation of program 
to text notation# Some symbols have different meaning in different portions 
of the text; however, their use in any particular instance should be clear 
from the text. 

Symbols used in text: 

b Length of wave crest between orthogonals 

C 1) Wave velocity (Phase velocity) 

2) In a Fresnel Integral, the real portion 

C^ Wave group velocity 

D Shoaling Coefficient 

d Depth of water - still level to bottom 

f Function of one or more variables, e.g. f(X,Y) 

2 

g 1) Acceleration of gravity (32.2 ft#/sec. ) 

2) Function of one or more variables, e.g. g(X,Y) 

H Wave height 

K, Refraction coefficient 

d 



Numbers in brackets refer to References listed in Section VIII. 
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K f Diffraction Coefficient 

K 2 tt/L 

L Wave length 

r Distance from end of breakwater to point (X, Y) 

S 1) Increment distance used when tracing a wave ray 

2) Imaginary portion of a Fresnel integral 
T Wave period 

t Time 

X Horizontal coordinate 

Y Vertical coordinate 

A (Delta) change 

O (Delta) change 

r) (eta) surface elevation 

7T pi - 3.14159 

0 Angular displacement 

-q Subscript "q" refers to deep water value of subscripted 

property 
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I c INTRODUCTION 



Waves and their effects constitute a major part of civil engineering. 
Countless books and papers have been written on wave generation, wave 
theory, wave propagation, and wave effects on coastal areas and coastal 
facilities. Waves, and especially ocean waves, are already complex in a 
deep water environment; once the waves enter coastal areas where interac- 
tion takes place with a bottom surface and with coastal facilities such as 
breakwaters, the problems become even more complex. 

In shallow water, wave speed is a function of water depth. As waves 
approach a shoreline, wave crests are bent through the process of refraction 
to conform to bottom contours. This refraction process causes changes in 
both wave height and direction of travel, both of which are essential fac- 
tors in any coastal engineering problem. 

Through the process of diffraction, waves and energy are propagated 
into geometric shadow areas behind breakwaters while reflected waves and 
energy move into regions beyond the ends of breakwaters. Wave patterns in 
diffraction zones are a complex of exponentially decreasing amplitude waves 
spread in an arc in the lee of the breakwater, while sets of waves alter- 
nately reinforce and cancel each other in an irregular pattern beyond the 
breakwater tip. 

From the second half of the last century through World War II, these 
problems were solved by empirical methods and through model studies. How- 
ever, due to military necessity for the conduct of operations across beaches, 
extensive research and analysis was conducted during the War, and resulted 
in major advances which formed the basis of the science as it now stands 
today . 
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At present, construction of wave refraction diagrams for use in 
engineering analysis requires manual drafting techniques such as those pre- 
sented by Johnson, O’Brien, and Isaacs [2], However, these methods are 
slow, tedious, and in many ways subjective. Although there have been 
various mathematical solutions for special cases of sloping, circular, 
shelf, and other analytical beaches, these solutions are of limited use to 
engineers in solving problems of natural beaches . 

Penney and Price [3] showed that the diffraction of light is also a 
solution to the water wave diffraction problem. Using Penney and Price’s 
theory as a basis, various diffraction diagrams have been published showing 
wave fronts and coefficients of diffraction for standard breakwater prob- 
lems. Due to the extreme complexity of the mathematics of the problem, 
various approximations are usually made to simplify the solutions. Results 
sometimes reflect these simplifications. For engineering studies, diffrac- 
tion analysis is performed through the use of overlays and tracings based 
on dimensionless plots which have been published. These plots are avail- 
able for standard diffraction problems at semi-infinite breakwaters and at 
breakwater gaps of various widths. 

In the interest of good engineering analysis of problems, and in view 
of the capabilities which are now available to the engineer through the 
use of electronic computers, a numerical analysis of these problems which 
is only limited by the basic theory is possible. The basic intent of 
this study was to develop a program system which would accept as input 
the same raw data which is available to an engineer; recognize any limita- 
tions which the engineer wished to impose on the problem; then carry out 
whatever sequence of commands the engineer wishes to use in the thorough 
analysis of any coastal area. 
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The problem of adapting the refraction process to a computer was 
first approached by Lt. G. M. Griswold of the U. S. Navy Weather Research 
Facility [4] in 1963, His efforts and later attempts by others have had 
varying degrees of success. Early programs often gave erroneous results; 
later programs worked better but did not calculate all the information 
needed by engineers; none of the programs were readily suited for use by 
practicing engineers. Analysis of the refraction problem therefore meant 
picking up where prior studies had left off. 

The major problem involved in refraction analysis is that of simula- 
ting bottom topography in a computer. This was the principle difficulty 
experienced with programs which were developed and tested by Harrison [5] . 

In view of the success which had been achieved by Snyder [6] at 
approximating planar curves through a method of continuous parabolic inter- 
polation, it was felt that the same principle, adapted to three-dimensional 
surfaces, offered an excellent yet simple means of surface fitting for 
refraction studies. The method would assure continuous quadratic surfaces 
at any point within the data structure. A method was therefore developed 
for use of three-dimensional continuous parabolic interpolation in tracing 
wave rays over a field of varying depths. The method was tested on slop- 
ing planar analytical beaches and gave results which agreed within 1% of 
the actual theoretical values. When tested on a natural beach, it gave 
results which very closely approximate the graphical procedures now being 
used. Although more extensive testing is desired, the system developed 
appears to offer a means of tracing wave rays which is readily suited to 
use in engineering analysis. The system is of as good or better accuracy 
than present methods, and yet can offer electronic computer capability for 
rapid economical solutions. 
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As best as could be determined, the problem of adapting the diffraction 
process to a computer has not been approached* Various methods have been 
developed for determining coefficients of diffraction using the original 
Penney and Price theories* Without too much difficulty, these methods 
are adaptable to numerical analysis „ For accurate engineering analysis 
however, it is also necessary to know the direction of wave travel at any 
given point in a diffraction zone* Although methods are available for cal- 
culating wave phases at any point in a diffraction zone, the only method 
at present to calculate wave direction of travel is through a graphical 
analysis of many points in an area surrounding the point in question. No 
means of directly calculating the direction of wave travel at a point in a 
diffraction zone has been published. 

Solutions of the diffraction problem therefore involved adapting 
existing diffraction methods to a numerical procedure for the computation 
of diffraction coefficients; and the algebraic analysis of a diffracted 
wave front to determine its direction of travel then adapting this analy- 
sis to a numerical procedure. Results obtained from the analysis and 
programs are comparable to those achieved through manual and graphical 
techniques . 

The refraction and diffraction programs which were developed were 
then incorporated along with an algebraic program to locate shorelines in- 
to a system which could readily be used by a practicing engineer. A 
problem-oriented language was used for operation command definition; input 
is the same basic data available to the engineer in the field; output is 
the same product information desired by an engineer in any analysis. 
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Use of this type of program could mean that, given necessary input 
data, engineering analysis of an area could be completed by one engineer 
within a matter of minutes rather than in a matter of months as is the 
case using older methods. 
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II. REFRACTION ANALYSIS 



A. Retraction Principles 

The velocity of a wave, C, varies with water depth in accordance with 
the basic relationship presented by Eagleson and Dean [7] and others 



gT 



2tt 



tanh 



27Td 

L 



( 1 ) 



where g is the acceleration of gravity, L is the wave length and d the 
water depth. With a constant period, velocity and wave length decrease 
as depth decreases. In deep water, 



and 



tanh ^=1,00 

L 



gT 

271- 



Looking at the crest of a long wave moving at an angle to the shore, 
it can be seen that in accordance with equation (1) the deep water portions 
of the crest move faster than the shallow water portions. This velocity 
variation causes the wave length to change and the wave to bend toward 
alignment with the contours. Refraction is the process whereby the di- 
rection of motion changes due to a change in wave velocity. 

Refraction results in changes in wave height, direction of travel, 
and change all characteristics of the wave except for its period which 
still remains constant. The extent of these changes depends on the 
bottom topography, the wave period, and original deep water direction of 
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travel. In performing refraction analysis, it is assumed that as a wave 
travels shoreward, no energy flows laterally along a crest; the energy trans- 
mitted between orthogonals to the wave crests will remain constant. The 
validity of this assumption will be discussed later. Using energy con- 
siderations it can readily be shown, Wiegel [8] (pp. 156), Beach Erosion 
Board [9], and Johnson [10], that 




where H is the wave height. C is the group velocity, and b is the 

S 

orthogonal spacing. The shoaling coefficient , D, is defined by 



D - 



which can be evaluated through the use of : 




C = n C 
g 



„ - i [1 + 4T,d/L 1 

sinh 47Td/L 



n = — in deep water 



This shoaling coefficient gives the ratio of the wave height in 
water at some shallower depth to the height of the wave in deep water 
with no refraction effects. 
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This study however will be more concerned with the derivation of the 



coefficient of refraction , K^, defined by: 




( 2 ) 



The coefficient of refraction gives the effect of wave crest bending, 
or changes in spacing of adjacent orthogonals, on the ratio of local wave 
height to the deep water wave height. 

The basic principle behind all refraction diagrams is Snell’s Law: 

sin A1 _ Cl 
sin A2 C2 

where A is the angle between a wave front and its respective bottom 
contour and C is as given by (1), This law states that, where bottom con- 
tours are parallel, the sine of the angle between the wave crest and the 
bottom contour is proportional to the velocity of the wave. When all con- 
tours are straight and parallel it makes no difference whether the change 
of depth is a continuous slope or a step, the change in wave length is 
only controlled by the deep water and shallow water wave lengths. 

Although equations have been derived for the solution of circular and 
parabolic beaches by Wiegel [8] (p. 160) and by Pocinki [11], for natural 
beaches, final wave direction is very dependent on intermediate contours 
and incremental procedures must be used for solutions. 

B, Present Methods of Refraction Analysis 



There are two methods in common use today; both are graphical, 
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The Wave Crest Method uses a graduated scale to plot wave advance from 
point to point along a crest, Starting in deep water, successive wave 
crests are plotted until the beach or end of the study area is reached* 

The Crestless Method plots wave orthogonals directly by determining 
shoreward deflections as orthogonals cross successive bottom contours, 

The amount of deflection is obtained from formulae derived from Snell 1 s Law, 
Both of these methods are described in considerable detail in various 
publications including Johnson, et al [2], Bretschneider [12], Dunham [13], 
and Beach Erosion Board [9]« In addition, an apparent improvement to the 
crestless method was developed by Arthur, Munk, and Isaacs [14] which gave 
closer to theoretical results when tested on concentric circular contours. 
Each method has its advantages and disadvantages and the method used in any 
particular situation depends on wave and topographic characteristics, and 
also the engineer /draf tsman ' s capabilities. Both methods, if applied and 
used properly, give results which are in reasonable agreement with aerial 
photographs of study areas as shown by Dunham [13]. 

However, there are two very basic problems with both systems. First 
of all, they depend very much on the operator’s skill, ability, and 
judgment. According to Johnson, et al [2], ’’Since bottom features which 
are comparatively small in respect. ito. the wave length do notr affect the 
wave appreciably” standard operating procedure is to ’’use judgement” and 
’’smooth out” irregularities. It is left to the engineer to define 
’’comparatively small” and ’’affect . „ .appreciably”. 

As shown by Dunham [13] and Pierson [15], different operators can get 
different results using the same system, the same operator can get marked- 
ly differing results using different systems. In particular, the wave 
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crest method may smooth out an area of orthogonal convergence and hide 
potentially dangerous situations, (Pierson [15]) „ 

A second basic problem is that a thorough study of a coastal area is 
extremely slow, time consuming, and thereby inherently expensive. Given 
that very small changes in direction or of period can cause extreme 
changes at a shoreline, it may be necessary to study small increments of 
wave periods for from 2 to 20 or more seconds over approximately 180 degree 
arc ranges, for many areas, at various tide stages . To do a reasonably 
complete job of studying an area can easily require 2 to 3 man months of 
effort; (Pierson, et al [16]). For complex study areas this time can in- 
crease considerably. For smaller low cost facilities or for those in which 
time is of critical importance, thorough engineering analysis quickly 
becomes impossiblec 

C. Previous Numerical Methods Used in Refraction Analysis 

In an effort to overcome the limitations and problems of graphical 
wave-ray diagrams, Lieutenant Gale M. Griswald [4], Then of the U, S. Navy 
Weather Research Facility, attempted in 1963 to develop a wave-ray trac- 
ing program based on the numerical methods originally proposed by Munk 
and Arthur [17]. His work, together with that of his associates, resulted 
in programs being written for the CDC 1604, the IBM 7090, and the Bendix 
G-15D computers, (Griswald [4] and [19]). However, these programs gave 
erroneous results in some circumstances apparently because of the bottom 
curve fitting procedures which were used. Results of these programs are 
shown by Harrison [5]. 

This work was then continued by the Army Coastal Engineering 
Research Center. Under contract to CERC, Dr« Wyman Harrison, of the 
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Virginia Institute of Marine Science, developed a method of constructing 
wave rays which appears operational (Harrison [5]) c However, for engin- 
eering use this program has many limitations. Because of the linear 
bottom curve fitting techniques used, it is not capable of calculating a 
coefficient of refraction as a ray moves shoreward. Independent passes 
of various programs are required to create input needed by succeeding 
programs, thereby resulting in rather complex operational procedures. These 
past numerical analyses and attempts to perform refraction analysis on 
computers were milestones. However, they did not go far enough towards 
giving the engineer all the information he needed in a simple efficient 
manner which he would be willing to use in practice. 

The basic approaches and theories behind these programs were sound 
however and optimum use of their methods and experience formed the basis 
for this study, 

D, Theoretical Considerations for Numerical Refraction 
a) Ray Curvature Equations 




Figure 1 
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As shown in Figure 1, the equations for a ray, S, making an angle, 



A, with the X axis are 



dX 

dt 



dS 

dt 



cos (A) 



dY 

dt 



dS 

dt 



sm (A) 



dS 

Letting — = C, we have 



dX 

dt 



C cos (A) 



(3) 



dY 

dt 



C sin (A) 



(4) 



Munk and Arthur [17] apply Fermat's principle'*’ to water waves and 
show that the curvature of the ray depends on the gradient of the velocity 
field normal to the ray'. Thus 

dA _ _ dC 

dt dn 

„ , . dS „ , . dA dA . dS 

Substituting = C and using — = — — we get 



dA = _ 1 dC 

dS C dn (5) 



Fermat's Principle states that light travels between two points along 
that path for which the travel time is a minimum. 
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dC 

We can now express the gradient — in a more convenient form by using 



the chain rule 



d _ _3__ 3x. 3__ _3y_ 

dn 3x 3n 3y 3n 



Noting from Figure 1 that we can show 



3x 

3n 



sin A 



3y 

“ - cos A 
3n 



and 



3_ 

3n 



3 3 

— (-sin A) + — (cos A) 
dx dy 



Defining a new curvature of the ray, , 



equal to 3A/3S we have from 



(5) and (6) 



K 

s 



3A 

3S 



1 

C 



, 3C . 3C n 

( sin — - cos A — ) 



(7) 



It is possible to express wave velocity and its derivatives from 
bottom depth and its derivatives. 



9C _ 9d/3x 
9x ” 9d/9C 



9C = 9d/9y 
3y 9d/9C 



( 8 ) 



(9) 



One cal calculate 9d/9C from equation (1) as follows: 
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tanh 





Setting 



K* = T/4n 




= 2CK' * i [In (1 + CK") -in ( 1-CK") ] 

44 - = K' [In (1 i- CK") -in (1-CK")] + CK 



[_£ 

1 1+CK" 1-CK" 



K^ ^ 



dC = K ' [ l- f cK"V 2 +ln (1+CK,,) _ln d-CK")] 



( 10 ) 



Given a means of evaluating bottom surface derivatives, equations 
(3) through (10) can be evaluated to find the curvature of a wave ortho- 
gonal at any point and the problem becomes one of tracing this ray of 
varying curvature from one point to another through a suitable program. 

b) Wave Intensity equations 

Defining a wave intensity factor, or more accurately a ray separation 
factor, (3, (Beta) by 



where b is the distance between adjacent rays, Munk and Arthur [19] de- 
rive the following equation for wave intensity along a refracted ray. 




(ID 




2 



( 12 ) 
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( 14 ) 



P(S) 

q (s) 



A r 1 3Ci A r 1 3C n / 10 \ 

= “ COS A [— - sin A [— (13) 



= sin 






2 

|^§]-2 sin A cos A [•£ 



3 2 C i 
3X3Y J 



2 A r l 
+ cos A [— 



3 2 C-; 



3 Y^ 



Factors dd/dC, 9C/dX, and 9C/9Y have previsously been defined by equations 
(10), (8), and (9), To evaluate 9^C/9X^, 9^C/9x9y, and 9^C/9Y^ we proceed 
as follows* 



9C dC , 9d 

9X V dd' l 9X ; 



3 2 C ,d 2 C . 3d. ,3dv £c . 3 2 d 

3 ^ = ( dI 7 3 l )( 3 ~ } dd “2 



3 2 C = _ d 2 d/dC 2 .3d 2 3 2 d/3X 2 

3X 7 '(dd/dC) 3 ^3x' + dd/dC 



Similarly 

2 2 2 2 ' 2 

9 C = _ d d/dC 3d 2 9 d/9Y 

W 7 (dd/dC) 3 l 3Y ; * dd/dC 



and 

3 2 C = _ d 2 d/dC 2 . 3d 3d 3 2 d/3x3Y 

3X3Y “ (dd/dC) 3 3Y 3X 3d/dC 



2 2 

To evaluate d d/dC we continue the differentiation of equation 

( 10 ). 



ft- K ' M-MltCK") - In (1-CK") ] +CK' [f^Tr - ^ 
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( 15 ) 




^ *'1 + CK") 2 + (1-CK") 2) ] 





K ’ ^l-CK" + 1-CK" 



K" . K" 



+ Ct (l-CK") 2 " (li-CK") 2 ^ 



dC 



Given a means of computing bottom surface derivatives, at any point, 
equations (13) through (15) can be evaluated. A means of solving (12) 
for 8 is then needed*, This is discussed m the succeeding section. 

Once 8 is found, the coefficient of refraction is obtained from 
equations (2) and (11) to give 



E. Program Considerations 
a) Basic Program Approach 

Solution of refraction problems involves finding some means of 
tracing a wave ray of known characteristics from deep through shoaling 
water. As explained earlier, at any point on a natural beach, the wave 
characteristics are very dependent on the bottom topography between the 
point and deep water, Procedures which allow for these intermediate 
effects, must be used for solutions * 

Incremental procedures used and tested by Griswald [4] and [19] and 
by Harrison [5] were basically sound. Therefore, similar techniques were 
used in the numerical process developed herein. 



K d = B- 2 



(16) 
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Using the known wave ray direction of travel and bottom characteristics 
at a given point, the ray curvature can be calculated using equation (7) . 
Using the initial direction of travel, and knowing the curvature, an es- 
timated location of the ray some incremental distance ahead can be 
determined* At this estimated location, it is possible to calculate a 
new ray curvature which can then be averaged with the curvature at the 
initial point 0 Using this average curvature, one can then revise the es- 
timated new location of the ray. At this newly estimated location, the 
curvature can again be calculated, a new average curvature determined, and 
the process can continue, until the change in curvature between successive 
position estimates reduces to within an acceptable limit. 
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It is possible to calculate all desired wave ray information at this 
new point including the direction of travel. From this new point, the 
next increment of advance can be determined using the same iterative pro- 
cedures. The process can continue until the wave ray runs up onto a 
beach, strikes a breakwater zone, or goes beyond the limits of the array e 

The program developed includes procedures which automatically check 
for shorelines, array limits, and intersection with a breakwater or a 
breakwater diffraction zone. 

To save unnecessary calculations, it is desirable to have the wave 
ray advance rapidly, using long intervals in deep water where bottom 
effects are small, and slowly, using shorter intervals, in shallow water 
where bottom effects are more important * The ray trace process therefore 
includes a scheme whereby, if a long interval is specified for deep water 
calculations, the interval will be reduced in successive stages as the 
ray curvature increases. 

These processes and the manner in which they are used in the program 
are explained in detail in Section IV and Appendix B. 

b) Depth matrix vs » velocity matrix 

Prior numerical methods for making wave-ray computations have been 
based on using an initial computer program to change an original depth 
array into a wave velocity array; then, as necessary, interpolating and 
deriving velocity values and derivatives from this new velocity matrix 
for direct input into the velocity/ordinate equations, such as (7), (13), 
and (14) . Making this change to a velocity matrix decreases the amount 
of calculation needed during the processes of tracing a ray and calculating 
refraction coefficients. It does this by eliminating the use of chain-rule 
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multiplication to convert velocity/depth and depth/ordinate derivatives 
into the necessary velocity/ordinate derivatives. 

However, a depth matrix is used in these programs for the following 
reasons : 

(a) An engineer could use mixes of wave periods during the analysis 
of an area without going through the lengthy procedure of resubmitting 
his original depth matrix and converting it to a new velocity matrix which 
corresponds to the changed wave period. 

(b) Water levels can be readily changed to simulate tidal 
fluctuations « 

(c) During possible expansion of the programs for further use in 
coastal analysis beyond refraction/diffraction studies, a depth matrix 
should have fewer limitations than a velocity matrix. 

c) Interpolation surfaces 

Various schemes for representing the bottom topography were tested 
during past analysis of the refraction problem by Harrison and Wilson [5]. 

A "f orced-cubic" interpolation scheme was originally used by 
Gr iswold-Mehr which derived a cubic surface of best fit to the velocity 
values at 12 grid intersections. This surface was forced to pass through 
the 4 values closest to the point in question and was the cubic of best 
fit for the next eight closest eight grid intersections. As a cubic, it 
permitted taking second derivatives for use in wave intensity calcula- 
tions. When this scheme was tested by Harrison and Wilson [5] it was 
found to give completely erroneous results in some cases of a wave ray 
approaching a beach at a shallow angle* 
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A "quadratic-interpolation" surface which was also tested by 
Harrison derived a quadratic surface of best fit to 12 grid intersections. 
This scheme gave good results except when tested within 2 grid units of 
the shoreline. Here it was found to give excessive curvatures, 

A "linear- interpolation" surface was then developed by Harrison and 
Wilson, and was based on the assumption that velocity values in a grid 
cell could be represented by a plane, This scheme gave good results up 
to within one grid unit of the shoreline but, as a linear surface, did 
not allow calculation of the derivatives which are needed for wave 
intensity calculations , 

For use by engineers, it was felt that calculation of wave intensity 
should be a minimum requirement of any program; this necessitates use of 
at least a quadratic surface, In view of the recent experience of others, 
a quadratic surface of best fit using least squares and possibly with 
heavily weighted central grid values appeared to be a logical method. 
However, this surface would not be continuous as the ray moved from one 
grid to another, Munk and Arthur [17] point out that for successful 
solution of equation (5) a continuous surface is necessary. 

It was found that excellent results at deriving continuous planar 
curves using parabolic interpolation had been achieved by Synder [6] in 
problems not related to refraction. 

Based on the geometrical device of linearly transforming the equa- 
tion of a parabola across an interpolation interval, a smooth curve can 
be developed which is continuous with no abrupt changes in slope even at 
data points. The function and its derivatives at the end of one interval 
are the same as at the start of the next. 
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The method was developed as a compromise between simple interpolation 
with a large input requirement and a more complex interpolation scheme 
with longer computation time £ When tested on long segments of a sine 
wave the method gave results which were almost as good as a cubic equationo 
Continuous parabolic interpolation appears to offer a curve fitting scheme 
that is little more difficult than for a simple parabola yet is actually 
a "third degree 11 equation* Its principal advantage however is the con- 
tinuous nature of the function which results. 




Continuous Parabolic Interpolation 
Figure 2 

Continuous parabolic interpolation works as shown in Fig. 2. Given 

4 data points, 1, 2, 3 and 4, the problem is one of finding the value at 

abcissa 0 of a smooth continuous curve which passes through all 4 data 

points. In order to calculate the value and slope at abcissa 0, a para- 

2 

bola of the form ax + bx + c is drawn through points 1, 2, and 3 and its 
projected value at the abcissa of point 4 is determined* An imaginary 
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point, 4 f , is located by linearly interpolating the difference between the 
projection of parabola 1-2-3 and the actual point 4 value proportional to 
the advance of point 0 from point 2 to point 3, A parabola through the 
imaginary point, 4 f , is used for calculating the curve value and its 
derivatives at point 0* This interpolated parabola is only valid at 
point 0, It can readily be seen however that as point 0 moves from point 
2 to point 3 a smooth transition is made from parabola 1-2-3 to parabola 
2-3-4 and a smooth trace of intermediate data values and derivatives is 
generated c 

Although only used in the past for curve fitting, it was felt that 
through the use of multiple parallel and orthogonal parabolas the same 
basic scheme could be adapted to definition of a 3-dimensional surface. 
Just as the scheme generated smooth curves in one plane it should generate 
smooth surfaces in 3 dimensions. 

The simplicity of the interpolating parabola method, especially 
when performed on a computer, and the fact that it met the other refrac- 
tion criteria of providing smooth curves and continuous derivatives at 
all points made it a best choice for a first try at surface fitting. 




Continuous Parabolic Interpolation 
for Surface Fitting 

Figure 3 
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For surface fitting, continuous parabolic interpolation was used as 
shown in Fig e 3* Using the parabolic interpolation techniques just de- 
scribed, 4 independent parabolas were used along vertical planes at A, B, 

? 

3d d d 

C, and D to calculate d, and — r for each of points (1, Y) , (2, Y) , 

Y 9Y^ 



(3, Y) , and (4, Y) a A parabola was then interpolated along plane E to 

determine (X, Y) values,, One "E" parabola based on 4 depth values gave 
2 

3d ad 

d, ’rrr, and — r* at (X, Y) * Another interpolated parabola based on the 
3d 3X 3d 3 2 d 

^ values gave an interplated 7 ^ 7 , and ^ or Y) and, similarly, one 

3 2 d s^d 

based on the 4 — — values gave - — ^ for (X, Y) . 

3Y 3Y 

Examination reveals that identical results would be obtained by 
reversing the X and Y order and first performing 4 interpolations along 
parallel X planes and the final interpolations along a Y plane. To con- 
firm this fact orthogonal tests were run on real data samples; as expected, 
both tests yielded identical results for (X, Y) and all its derivatives* 

By geometry, it is readily possible to reduce the interpolation 
procedures described above to simple functions. Using the notation of 
Fig. 2, but with the 4 equally spaced data values designated Al, A2 , A3, 
and A4; and with D defined as 6 /A, which is always measured from A2 
towards A3; one can solve for the value of AO. 



AO = A2 + (D/2) (A3 - Al) - (D 2 /2) (A4 - 4A3 + 5A2 - 2A1) 
+ (D 3 / 2) (A4 - 3A3 + 3A2 - Al) 



Its first and second derivatives respectively are 

AO' = (A3 - Al)/2 - (D/2) (A4 - 5A3 + 7A2 - 3A1) 
+ (D 2 ) (A4 - 3 A3 + 3A2 - Al) 

AO' ' = A3 - 2A2 + Al + (D) (A4 - 3A3 + 3A2 - Al) 
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These equations used as functions in a computer program greatly 
facilitate the solution of continuous parabolic interpolation problems* 



d) Solution of the Wave Intensity formula 

Various methods for solving the wave intensity formula equation (12) 



D 2 f3 , n DB , n R _ n 

D s 2 P DS q 6 



are discussed by Munk and Arthur [17] • These methods include: 

(a) Using average constant values for p and q within each interval* 
The equation can be solved and the incremental solutions joined in a con- 
tinuous curve* This solution appeared complex and not readily suited for 
incorporation into the other refraction programs, 

(b) A "WKB M method involves transformation to forms used in quantum 
mechanics. This method appeared more complex than (a) preceeding and 
was ruled out in view of the simpler methods which were available, 

(c) Analogue computer methods could be used to solve this problem 
but were ruled out as beyond the scope of the problem* 

(d) Kelvin 1 s method involves approximating the integral curve by 
fitting together circular arcs. The method is readily adaptable to the 
associated refraction programs <, This method was tested by Griswold [4] 

to calculate wave intensity for an analytic field of wave velocity values. 
Wave intensity was then compared with actual theoretical values. The 
greatest deviation from theoretical was 1.6 per cent. 

In addition, Griswold [4] used a simple finite difference method. 
Using 3 successive values of designated BN, Bl, and B2, at equal 
spacing, S, along a curve, one can find by geometry 
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( 17 ) 



£ (BN - 2B1 + B2) / (S 2 ) 
dS 

4| S (B2 - BN) / 2S (18) 



Substitution of (17) and (18) into (12) and solution of the resulting 
equation gives 

B2 ■ t V+ffi B1 * BN < l9 > 

Application of this formula requires that, for three equally spaces 
points, we know Bl, p, and q at the center point, and that we know BN at 
the initial point 6 The formula calculates B2 at the third point. This 
method requires no iterations, and simply involves projecting one Beta 
value ahead to the next succeeding point, B2 0 After moving down the ray 
trace, B2 is then used for the next Bl value and Bl becomes BN. This 
method was also tested by Griswold [4] against the same analytic field 
used for the Kelvin method, and was found to give errors of less than 
2.2 per cent. 

In view of the much greater simplicity of the finite difference 

method and the fact that it appeared to give results almost as good as 

this method 

the far more complex, Kelvin’s method/was used in the refraction program. 

F . Tests applied to Ray Trace Program 

The theory of the Ray Trace program was tested using two methods. 
First, it was tested on analytical beaches of uniform slope from 
deep water up and onto a beach* These tests could be readily verified by 
means of Snell’s Law: 
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sin Al sin A2 



( 20 ) 



Cl C2 



and a Beta expression derived from Snell’s Law which is given by Munk and 
Arthur [17] 



sm A2 
sin Al 



( 21 ) 



One of these tests, which could be considered a typical case is shown 
in Appendix D, This test shows a wave ray on a plane beach varying in 
depth from 60 feet to zero. The ray has an initial angle of 45 degrees 
and a period of 4 seconds* Its error analysis, based on equations (20) 
and (21) is shown in Table 1. This and similar tests on plane beaches 
gave results which agree within 1 percent of the theoretical values for 
both direction of travel and Beta value. 
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Secondly, the method was tested on a section of natural beach at 
Virginia Beach, Virginia. This test was selected due to the fact that 
previous refraction methods were tested on this beach by Harrison and 
Wilson [5], and results achieved could therefore be tested both against 
graphically constructed wave rays and against other computer programs. 

The most accurate method developed by Harrison and Wilson [5] for 
this beach involved using interpolation with a linear surface of best fit. 
Therefore, the results of the continuous parabolic interpolation tests are 
compared to the linear surface method and the more common graphically con- 
structed rays. The comparison is shown on Fig. 4. 

The accuracy of these comparisons was limited by the fact that the 
input data used by Harrison was not available. It was necessary to es- 
timate input data and his solutions from a 3" by 4" drawing. This could 
have lead to some differences in results but it was felt that data were 
suitable for initial testing. 

These limited results make it appear that the method developed is 
comparable to the graphical scheme and linear surface method. 

It is also important to remember that none of the methods can be 
verified under natural conditions for reasons discussed elsewhere in this 
report. Although graphical methods are most commonly used in practice 
today, different graphical methods would give different results, different 
draftsmen using the same method would give different results. 

One cannot say with certainty which of the three methods most closely 
approximates what would actually happen at this beach. 

In addition, many specialized minor tests were given the refraction 
program to assure that special features such as the breakwater intersec- 
tion check and the gradual reduction of the increment interval do in fact 
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Refraction Test - Natural Beach 



work as programmed. These tests gave what amount to Yes/No results,, 

These tests are not enumerated here in view of their lack of significance 
once they work as desired. 

Go Program Limitations 

a. Limitations of the Numerical Methods 

This program is limited in various ways as a result of the numerical 
methods used a Wave ray tracing is based on tracing a wave ray using a 
series of finite intervals. At any point on a natural beach, wave charac- 
teristics are dependent on bottom topography between that point and deep 
water, Any finite interval techniques used may skip over intermittent 
bottom effects. Extensive testing is needed to determine the effects of 
this interval technique* When various intervals were tested on the pro- 
grams developed, the effects were very small; but more testing is needed 
especially on complex bottoms. 

This program is still limited by the ability of parabolic interpola- 
tion to accurately represent the bottom surface. In cases of extremely 
complex bottom topography, it is necessary that grid sizes be small enough 
that bottom features can be reasonably well represented by a quadratic 
surface over a span of 2 grid units. This means that over relatively 
smooth surfaces, large grids are possible; over an irregular bottom, 
smaller grids should be used, 

The solution of equation (12) is obtained by a finite difference 
scheme. Although this scheme has been tested against analytical beaches 
by Griswold [4] , it should be tested against natural beaches providing 
some suitable means of verifying results can be obtained* 
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b, Theoretical Limitations 



Most of the same original basic assumptions which were involved in the 
graphical construction of refraction diagrams are still present in the nu- 
merical methods. No new theory has been added; the final results are still 
limited by the accuracy of this theory c 

Refraction Theory is based on the assumption that waves are long 
crested* Short crested waves, especially if they cross contours at a high 
angle of incidence, will have refraction coefficients which differ from 
theoretical, Wiegel (8, p 172) « 

Although no limits are set, beach slope, roughness, and reflections 
can effect refraction on natural beaches because of their effect on energy 
distribution assumptions used in refraction analysis e 

The assumption that no energy travels laterally along a crest does 
not hold when orthogonals bend sharply or form a caustic envelope. A cer- 
tain amount of energy is transmitted across orthogonals in cases of complex 
bottoms and where sizeable variations in wave height occur along a crest e 
Waves do not actually become infinitely high as indicated by a caustic 
envelope; sometimes they peak and break but usually only a chaotic appear- 
ance results according to Pierson, [15] a Similarly, crossed orthogonals 
may indicate a severence of a wave train and subsequent crossed wave 
crests, but this is an area which is still full of uncertainties. 

Much is still left to the engineer T s experience and judgement. 

There is no precise method for smoothing observed contours to that point 
at which they will still represent wave refraction but where neglect of 
minor bottom irregularies is justified. In all cases, accuracy of output 
is limited by accuracy of the depth data available for input. 
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Ill, DIFFRACTION 



Ac Principles of Diffraction 
a) General 

Diffraction is the process whereby wave energy is transferred 
laterally along a wave crest and thereby propagated into areas of geo- 
metric shadow behind Impervious barriers „ It is similar to the diffrac- 
tion of light that takes place when light is transmitted sideways behind 
an object so that the area behind the object is slightly lit c Sound and 
electromagmetic waves experience similar effects. 

Diffraction can best be explained by Huygen's Principle which 
states that each point of a wave front is a source of energy and that 
it radiates outward equally in all foreward directions (Russell and 
MacMillan [18]), The wave motion at any point is the sum of the 
motions induced by all energy sources behind it. In this manner, the 
"end" of a wave which has been interrupted by a breakwater will act as 
a source of energy and, in the lee of the breakwater, will spread out in 
a circular arc of exponentially decreasing amplitude, (Wiegel [8], 
p 0 100). That wave which is reflected off the face of the breakwater 
also has an "end" which is a source of energy for those regions beyond 
the breakwater tip f This results in two sets of waves, one advancing 
and one radiating, which alternately reinforce and cancel each other so 
as to cause very irregular wave heights in this region beyond the tip. 

As a wave approaches a gap, diffraction takes place in both direc- 
tions as it passes through the gap and two sets of radial waves result 
from the reflected waves 0 Addition of the advancing waves and two 
radial waves in the lee of the breakwater results in an extremely 
complicated pattern, 
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Penney and Price [3] show that the distribution, amplitude, and 
phase of surface waves on a sheet water of any uniform depth can be de- 
scribed by a complex wave function, F(X, Y) . Through the use of 
Sommerf eld f s solution for the diffraction of light, polarized in a plane 
parallel to the edge of a semi-infinite screen, they are able to resolve 
the water wave diffraction problem. 

The complete solution of the water wave diffraction problem has 
been given in considerable detail by Penney and Price [3] , and has since 
been summarized in various forms by many authors, including Putnam and 
Arthur [20], Weigel [21] and [ 8 ] ( p 180), and Bretschneider [12]. The 
reader is referred to these publications for the details and reasoning 
behind equations which are included in summary form only in this study. 



b) Semi-infinite breakwater 

Penney and Price [3] show that the equation for the surface 
elevation, n, in cylindrical coordinates (see Fig. 4) can be stated as 



n 



AikC 

g 



ikCt 

e 



cosh (kd) F(r,9) 



where A is a constant and k = 2 tt/L. F(r,0) is a function which must 
satisfy 



9 2 F 1 9F , 1 9 2 F 

2 r 3r + 2 2 

9r r ae 



+ k F = 0 



( 22 ) 



The diffraction coefficient, K f , is defined as the ratio of wave 
height in the area affected by diffraction to the incident wave height 
and is given by the modulus of F(r,0) or 

K' = | F(r ,0) | 
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In addition, the wave pattern or phase of the waves can be 
determined from the argument of F(r,0) which is denoted by arg F(r,0). 

The solution of the diffraction problem is therefore the solution 
of equation (22) , subject to the boundary condition that the normal 
component of the fluid velocity is zero along the breakwater „ 




Breakwater 



(S) 




Incident Wave Travel 



Notations for Semi-inf inite Breakwater 



Figure 4 



Using the notation of Fig. 4, and defining 



a 




(23) 




(24) 



Penney and Price show that in region "S" (i.e., 0 £ 0 < 0^) 



F (r , 0) = f (a) e 



-ikr cos (0-0 n ) 



Cr +f (a* ) e 



-ikr cos (0+0 A ) 
: 0 



(25) 
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and in region "Q" (i.e., 0^ < 0 < (0 q+tt) ) 



-lkr cos (0-0„) , , N ~ikr cos (0-0 A ) 

F(r,0) = e 0 -f(o) e 0 



+ f(o’) e" ikr COS C0+6 O ) 



( 26 ) 



where f(o) is defined by 



f (o) = \ (1 + i) J e 2 77111 du 



(27) 



Equation (27) can be evaluated through the use of Fresnel Integrals 
which are defined by 



C - iS = I e 

0 



■/ 



-ITT 11/2 



du 



or alternately, 



C(c ) 



-I 



1 2 

COS — IT u du 



(28) 



S (0 ) 



o 

f . 1 . 2 

-I sin — tt u 



du 



(29) 



Substitution into (27) gives 



f(o) = f [(1 + C + S) - i (S-C) ] 



Using the fact that 



f(-o) = 1 - f(a) 
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we can define 



f(~o) = U (a) + iW(a) (30) 

where 

U = J (1 - C - S) (31) 

W = y ( S - C) (32) 



Setting 



^ -ikr cos 

f = e 


0 

CD 

1 

CD 


t(o) 


(33) 


-ikr cos 

g = e 


(e+e o ) 


f (°') 


(34) 


f(-a) = U1 t 


i W1 




(35) 


f(-a T ) = U2 4 


i W2 




(36) 



it can be shown algebraically that 

f = U1 cos (kr cos (0-0^)) + U2 cos (kr cos (0+0^)) 

+ W1 sin (kr cos (0-0^)) + W2 sin (kr cos (0+0 q) ) 

g = i {W1 cos (kr cos (0-0 q)) 4- W2 cos (kr cos (0+0g) ) 

- U1 sin (kr cos (0-0^)) -U2 sin (kr cos (0+0^))} 



(37) 



(38) 
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Returning to our original equations it is seen from (25) that in 
region S 



F(r ,0) = f + g 
while in region Q, (26) gives 

F(r,6) - e‘ lkr cos - £ + g 

In either case, F(r,0) can be expressed as 
F(r ,6) = A + i B 



(39) 



(40) 



(41) 



hence 



K' = | F(r ,0) | = (A 2 +B 2 ) 1/2 (42) 

and the phase difference, <J>, can be expressed as 

<j> = tan' 1 (B/A) (43) 

Therefore by evaluating a and a f , solving for C and S in (28) 
and (29) , and making the transformations indicated by equations (31) 
through (43) , we are able to solve for the diffraction coefficient and 
phase difference at any point in the diffraction zone. 
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c) Breakwater Gap 



Using the following notation 




Breakwater Gap Notation 
Figure 5 



Penney and Price [3] also show that for waves entering a breakwater 
gap at right angles (0 q = 90^) the diffraction pattern is described by 



F(r , 0) = f x + % 1 ~ f 2 + S 2 


for 


region S 


(44a) 


F(r , 0) = -£ 1 +g x + f 2 + g 2 


for 


region T 


(44b) 


F(r ,0) = e lkx -f x + g x - f 2 + g 2 


for 


region Q 


(44c) 



where the subscript indicates the breakwater wing from which the func- 
tion is measured. 

This solution is only valid, however, for waves entering normal 
to the gap« Upon entering a gap at an oblique angle, a given wave 
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reaches one tip first and diffracts there before the remainder of the 
wave reaches the other tip where it will also diffract. Therefore the two 
wave patterns are out of phase unless the wave crest happens to reach 
the second tip an exact multiple of wave periods after reaching the 
first © Equations (44) do not, in general, give valid results with 
oblique waves. However, an approximate solution for oblique waves is 
possible through the use of techniques which will be discussed later 0 

B. Present Methods of Diffraction Analysis 

It is apparent that the solution of the diffraction equations is 
a lengthy procedure when done by manual methods. Fortunately, however, 
they can be solved independently of wave speed and wave period; and 
solved in terms of wave length, L« 

When plotted in terms of X/L and Y/L coordinates, diffraction 
charts become dimensionless plots of diffraction coefficients which 
can then be applied to any water wave diffraction problem in water of 
uniform depth regardless of wave period or wave speed. These solutions 
have been completed and the results have been included with virtually 
every diffraction publication, including Wiegel [8](p> 180), 

Bretschneider [12], and Beach Erosion Board [ 9 ] ® For the engineer in 
the field, it then becomes simply a matter of picking one of these 
dimensionless plots, expanding or reducing the scale so as to fit his 
particular situation, and with the use of a properly scaled overlay, 
analyze the particular diffraction situation at hand. Different over- 
lays can be made to represent each wave period to be studied. 

The same procedure applies for studying diffraction at breakwater 
gaps. In these cases, the ratio of gap width, b, to wave length con- 
stitutes a second independent parameter; hence, a series of dimensionless 



plots are required* These have been published by many authors, including 
Johnson [24], and Wiegel (3, p 189), for differing b/L ratios from 1 to 
5, Beyond 5 wave lengths it is felt that the 2 breakwater tips which 
form the gap are far enough apart that they can be considered independent 
for all practical purposes * 

Plots of breakwater gaps for waves approaching at an angle are 
available (Weigel (8, p 193) and Johnson [10]) based on the theories of 
Morse and Rubenstein, which were developed by Carr and Stelzride [25] * 

The calculations for these plots are based on series of Mathieu functions 
which give exact theoretical solutions for incident waves of any angle 
of approach* 

In addition, Johnson [24] showed that a modified gap width as shown 
in Fig. 6 could be substituted for the actual gap and results achieved 
which compare very closely to the exact theory of Morse and Rubenstein* 
This modified gap solution allows the use of Penney and Price equations 
which were developed earlier. 




Modified Breakwater Gap 
Wave Approaching at Oblique Angle 

Figure 6 
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In areas sufficiently removed from the breakwater, simplifying 
approximations can be made for some terms of the various formulae. 

These approximations have been verified experimentally by Putnam and 
Arthur [20] e However, in regions close to a breakwater a complete 
solution is needed. 

In addition to Putnam and Arthur [20] , diffraction theory has been 
verified experimentally by Blue and Johnson [26]. Their experiments 
showed that actual amplitudes agreed very closely with theory in the lee 
of breakwaters and were somewhat less, but within 10 percent of, theore- 
tical beyond the breakwater tip in the unsheltered zone. This results in 
generally conservative results, i.e. , overdesigned facilities, if the 
theoretical results are used. 

In addition to coefficients of diffraction, it is also necessary 
to know the direction of wave travel at any given point in a diffraction 
zone. From equation (43) , it is readily possible to calculate the phase 
difference between the diffracted wave and the wave front of the original 
incident wave. With this as a basis, it is possible to calculate a 
number of phase differences within a zone, then by connecting zones of 
equal phases determine actual wave crest positions and direction of 
travel. However, no published method has been found which would di- 
rectly give a direction of wave travel for a given point. 

C, Previous Numerical Methods 

No publications were found concerning numerical analysis of water 
wave diffraction problems. However, the present formulation of the 
diffraction problem is readily susceptible to calculation of coeffi- 
cients of diffraction in their present form, hence it is possible that 
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portions of the problem have been done by computer as an aid to graphical 
analysis but were not considered worthy of publication. 

D. Theoretical Considerations for Numerical Diffraction 

a) Coefficient of Diffraction, K f 

Except for evaluation of Fresnel Integrals, equations (23) through 
(42) can be solved on a computer in essentially their present form 0 In 
the programs developed, a completely random orientation of breakwaters 
and incident waves was desired. This and a desire for reasonable program 
efficiency and memory size cause substantial bookkeeping problems but 
create no need for any new methods or procedures. 

As shown in Appendix A, a series solution was developed for the 
evaluation of the Fresnal Integrals, This was used for sigma values 
less than 3.0. For larger values, sine/cosine approximations given by 
Sparrow [23] were used. If used with sigma values greater than approxi- 
mately 3.0, the developed series solution requires greater precision 
than normally used on most computers. Double precision could be used 
but tests indicated that satisfactory evaluation, good to within 0.1 
percent of the Integral, could be achieved by using the two methods., 

A third approximation formula is given by Jahnke and Ende [22] which 
is suitable for evaluation of values in the intermediate range between 
the series solution and the sine/cosine approximations but this was 
not felt necessary. 

b) Direction of Travel of Wave Crest 

As noted earlier in equation (43), the phase, <j>, of the diffracted 
wave is given by (J> = tan ^ (B/A) 
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ir< of Travel 



Incident 
Wave 




“V 



Wave Front Nomenclature 



Figure 7 



The equation of the wave crest is given by 
(f) (z,x) = Const 

The normal to the wave crest is a vector which is given by 

«»-■§£■? + §£ t 

in which “x and k were the unit vectors in the x and z directions 
respectively. Therefore, representing the slope of the normal as 



The partial derivatives in (44) can be expressed in terms of A and 



B. Using (43) this gives 



3<t> 1 2 f I 9B B_ 3A, 

3z 1+(B/A) A dz ~ 2 3z 

A 



li _ 1 r l 3B B_ 3A. 

3x 1+(B/A) A 3x 2 3x J 

A 



Hence (44) becomes 



[— J = 
l 3Z J _L 



A — - 

A 3X 



3B 

3Z 



- B 



3A 

3X 

3A 

3Z 



(45) 



(46) 



(47) 



Equation 47 is our desired solution for the direction of wave 
travel. 

Given A and B from equations (39) or (40) , we must now evaluate 
3A/3X, 3B/3X, 3A/3Z and 3B/3Z. Fox these we will need the derivatives 
3r/3X, 3r/3Z, 30/3X, and 30/3Z. For the breakwater along the +Z axis, 



Therefore , 



x 




x = r sin 0 

z =• r cos 0 

tan 0 = — 
z 



3r 

3x 



OL a n 

— = sin 0 



3r 

dz 



= cos 0 
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96 _ cos 9 
9x r 



96 _ _ sin 6 
9z r 



For the breakwater along the -z axis 




z = -r cos 6 

tan 6 = - — 
z 



Therefore , 



9r 

3x 



sin 6 



9r 
9 z 



-cos 6 



98 _ cos 6 
3x r 

96 _ sin 6 
3z r 



If we consider X = f(r,9), the chain rule gives 



9B 9B 9r 9B 99 

9x 9r 9x 96 9x 



which becomes 



9B 9B . „ . 9B cos 6 

T— = T— sin 6 + TT 

9x 9r 96 r 
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Similarly, one gets 



9A 

3x 



3A . 0 

sm 0 + 
dr 



3 A cos 9 
36 r 



Using the same techniques we can evaluate and but in this case, 

results differ with the breakwater orientation. If the breakwater is 
along the +z axis, 



3B 

3z 



3B 0 

7C — COS 0 

3r 



3B sin 9 
30 r 



3A 

3z 



3A 

COS 0 
3r 



3A sin 6 
30 r 



while if it is along the -z axis 



3B 3B „ . 3B sin 0 

— = - t: cos 0 + rr 

3z 3r 30 r 



3A 3A „ , 3A sin 0 

7 — = - cos 9 + Tq 

3z 3r 30 r 



Still using the notation of Fig. 7 the following convenient notations 
are used: 



k = 2tt/ CT 



FM = kr cos (0-0^) 

3FM , /A A \ 

A — = k cos (0-0 n ) 
3r 0 



3FM , . /a a \ 

30-=- kr sin ( 0 - e O ) 



FP = kr cos (0+0^) 



| FP = k cos (0+0 n ) 
3r 0 



IF “ - kr sin (e+6 0 ) 
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f 




Introducing this notation into equations (39) and (41), in region S, one 
has 



A = U1 cos(FM) + U2 cos(FP) + W1 sin(FM) + W2 sin(FP) 



(48) 



Taking derivatives 



|A , |!!i cos(FM) - HI sin(FM) eos(FP) -U2 .in(FP) 

dr dr dr dr dr 



+ sin(FM) + W1 cos(FM) 4- sin(FP) + W2 cos(FP) ~ 

3r 3r 3r 3r 



ff - cos(FM) - U1 sin(FM) + |]p- cos(FP) - U2 sin(FP) ||^ 



+ ||p sin(FM) + W1 cos(FM) ||^ + |^p- sin(FP) + W2 cos(FP) ||^ 



Similarly , 



B = W1 cos(FM) + W2 cos(FP) - U1 sin(FM) - U2 sin(FP) 



(49) 



3B 3W1 ,™„3FM , 3W2 f . f . 3FP 

■tt— = r cos(FM) - W1 sm(FM)-r h t cos(FP) - W2 sin(FP) r 

dr dr dr dr dr 



f^- sin(FM) - U1 cos(FM) - f^- sin(FP) - U2 cos(FP) f^- 

3r 3r 3r 3r 



3B 3W1 3FM . 3W2 ,_ x „„ . 3FP 

Jq = Tg— cos(FM) - W1 sin(FM) cos(FP) - W2 sm(FP) 



sin(FM) - U1 cos(FM) |—£ - sin(FP) - U2 cos(FP) ||^- 



Similar equations for A, B, and their derivatives apply in region 
Q or at the breakwater gap and may be derived from equations (40) and 
(41), or (44) and (41), respectively. These equations are not detailed 
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in this report in view of the fact that they offer nothing new in the 
way of knowledge and are simply exercises in bookkeeping. 

For evaluation of these last groups of derivatives, one must deter- 
mine Ul, U2, Wl, W2 and their derivatives. Consider first the evaluation 
of Ul and Wl. We recall from (31), (32), and (35) 

Ul = (1-S-C) 

Wl = J (S-C) 

Taking derivatives 

3U1 1 ,3S 3C. 

~ ' 2 ( 37 + 37 ; 

3W1 = 1_ .3S_ _ 3C, 

3 r 2 '3r 3r' 

3U1 1_ 3S. j. 9C, 

30 2 '36 Se’ 

3W1 _ 1 , 3S. 3C. 

36 2 ^36 36 ' 

To evaluate derivatives of S and C it is convenient to set 




which gives 



3_t 

3r 




( 50 ) 
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Equation (29) therefore gives 



s(t 1 ) 




sm c 




dS _ dS_ d_c _ d_t_ . sin t 
dr dt dr dr ^ 



substituting (50) 



dS 

dr 



2 

tth 

sin 



da 

dr 



(■^ keeps sign of o) 



Also from equation (28) we derive 



c(t x ) 




cos t 
s/ 2irt 



dt 



dC = dC dt = dt 
dr dt dr dr 



JTirt 



dC 

dr 



2 

✓ TO . 

cos (~ 2 ~) 



da 

dr 




keeps sign of 



o) 



Using the same techniques we can evaluate 0 derivatives of S and C and 
get 



dS 

d0 




do 

d0 



( lf kee P s 



sign of a) 



dC 

de 



r v 0 \ 
COS (— 7T") 



d.1 

d6 



kee P s 



sign of o) 
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From (23) we can derive 



o = 



3o 

3r 



2 JF sin 2 < 6 - 9 0 ) 

/f? Sln 2 ^ 6-6 0 ) 




3o 

7e 




cos y (e-e 0 ) 



We evaluate U2 , W2 and their derivatives similarly, using d* and its 
derivatives, i.e. 




Therefore, given values of r, 0, 0 q, and k, we can back substitute 
through the preceeding series of equations to ultimately evaluate the 
direction of wave travel as given by equation (47) . 



c) Alternate Methods Considered 

The two methods used in diffraction analysis by engineers have been 
the Penney-Price and Morse Rubenstein methods. When evaluated by Carr 
and Stelzride [25], both were found to have advantages and disadvan- 
tages but the Penney-Price method was found to be the better* Ever 
since, the Penney-Price has been used almost exclusively except in the 
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evaluation of breakwater gaps where the more exact theoretical solutions 
of Morse Rubenstein have often been used. Principally due to the fact 
that more testing had been accomplished and much more information was 
available in publications on the Penney-Price methods, this approach was 
used in the analysis 0 

Once committed to this method, when it was later decided to add 
solutions of the breakwater gap problem to the program package, the 
Penney-Price solutions using the imaginary breakwater gap developed by 
Johnson were used in the interests of simplicity* Although not an exact 
solution when applied in the vicinity of the breakwater wings, in the 
hands of an experienced engineer these areas should cause no real problems. 

Many simplifications which omit negligible terms in various re- 
gions of a diffraction zone could be made c Due to the unknown effect 
which these approximations would have in the determination of the direc- 
tion of wave travel, an exact solution was used throughout * The actual 
computer time saved by making approximations would be very small and was 
not felt to be worth the chance of having to rewrite programs from 
approximate to more exact solutions if the approximations did not work* 
Working from an exact solution however, approximations could be added 
at a later time and results readily verified against the original 
solution. 

E. Tests applied 

To test the accuracy of these programs, they were used to calcu- 
late diffraction coefficients and directions of wave travel at many 
random points using various breakwater coordinates and various direc- 
tions of travel for the incident wave* Results were then checked 
against graphical solutions which have been published. 
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Problems were encountered in this verification however, due to the 
fact that virtually every published refraction diagram was based on 
approximations of one kind or another (usually not mentioned in the 
write up) which had a substantial effect on the quality of the plot. In 
fact, no two published diffraction diagrams from different authors 
appeared to agreee with each other as a result of these various approxi- 
mations o The program results did not completely agree with any published 
diffraction diagram c However, in the calculation of the value for a 
coefficient of diffraction, good correlation was found between the program 
and the diffraction chart in Wiegel (8, p,183). Also, in the calculation 
of the direction of wave front travel, good correlation was found with 
the diffraction chart published by the Beach Erosion Board (9, p. 38) . 

Table 2 compares program results with diffraction coefficients inter- 
polated from Wiegel and with directions of wave travel interpolated 
from the Beach Erosion Board, It should be remembered that in both 
cases values interpolated off the charts are no better than the accuracy 
with which the chart was drawn or one’s ability to interpolate off 
them. In addition, analysis of the basic equations, especially in the 
region around and beyond the breakwater tip, where the most discrepan- 
cies are found between published charts, indicates that the results 
achieved are in fact reasonable and in good agreement with the complete 
theoretical solution. 

F. Limitations of Numerical Diffraction Programs 

In that these programs are built on the same theoretical basis as 
the past plots, they suffer the same limitations. No new theories or 
improved results can be claimed over a completely accurate plot. The 
program 1 s chief asset is its ability to produce results in a small 
fraction of the time required by graphical methods. 



- 60 - 



X 

CU 

£ 

CO 



> 

G 

J-i 

H 



P 

CL. 



mCMCNmCMOLDCM 

^^OOOOONONOOON 



(«::::: 

t> 

< 

<t <t <r 

4_J c . 0 

orrrrr m o cm 

^ o> cr> 



> 

CO 

M 





H 




00 


CM 


'vD 


VO 


i—i 








o 


r^. 


i—i 


m 


o 


o 


O 




u 




00 


co 


^D 






00 


I—l 




CO 


vD 


vD 


co 




o 


00 


o 


CM 


CO 


4-1 


o 


o 


0 


• 


o 


0 


o 


• 


0 


o 


c 


o 


o 


• 


« 


o 


• 


CU 


o 


co 


1—1 


co 


m 


CM 


1—1 




i— 1 


r- 






CO 


o 


in 




o 




H 






vD 


00 


00 






00 


a\ 


CO 


CO 


i—l 


CO 


<r 


<3- 




ON 


ON 



CM 

<U 

i—l 

rQ 

CO 

H 



co 

to 

i—l 

CO 

$ 

g 

o 



a 

(0 

J-i 



CU - 


i — 1 


vD 


CO 


CM 4 - 


r^co 






m 






co 


dD <T 


m<r 


ao ^ 


i—l 


CM 


m 


in rH 


O ON 


o 


o 


m 


CO 


o 


m 


O co 


00 00 


a; 


• 


C 


■ 


• • 


o o 


o 


tt 


© 


o 


e> 


o 


o • 


o • 




O 


o 


o 


O i—l 


r — 1 O 


rH 


1—1 


o 


o 


i — 1 


o 


o o 


o o 







I—l 


CJN 


o 


00 


o 






<3- 


oc 


00 


i—i 




CD 


ON 




00 




X 






4-» 


o 


vD 


CM 


CM 


o 


VsD 


CO 




o 


<T 


CM 


m 


co 


00 


CO 


<r 


00 




cu 






CO - 


1—1 


CM 


m 


m 


*— 1 


O 




o 


1—1 


m 


CM 


X\ 


m 


O 


CM 


00 


00 




■u 






cu & 


o 


• 


o 


• 


o 


• 


• 


• 


© 


O 


o 


o 








o 






o 






H 


o 


o 


o 


o 


r— 1 


i—l 


o 


i— 1 


1—1 


o 


O 


O 


o 


O 


O 


o 




G 












































J-I 










































• 


cu 










































XI 


£ 










































cu 


u 












































o 










































o 












































g 


CO 












































CO 






y— N 




































cu 


cu 






CO 




































CO 


1 — 1 






X 


CM 


-d- 


o 


o 


*<r 


-d- 


CO 


m 


<3- 


o 


•<r 


CO 


o 


O 


CM 


o 


cD 




G 






U vH 


CO 


CM 


o 


o 




<r 


00 


CM 


•<r 


o 


<3- 


00 


o 


O 


CO 


o 


O 


£ 


p 






J-I 


































o 


J-i 








ao 


r-^ 


CM 


•<r 


00 


<r 




CM 


00 


sr 


CO 


<3- 


m 


00 


00 


CD 


m 


CM 


CU 


00 






'w' 




































£ 


CO 










































4-» 












































o 


0 












































CL 










































CO 












































CO 


1—^ 








< 


D 
































cu 


ON 








o 


CO 


o 


o 


VD 


o 


O 


<r 


CD 


o 


co 


cD 


o 


O 


m 


o 


-d" 


1—1 










o 


>0- 


o 


o 


m 


o 


O 


o 


m 


o 




CTC 


o 


o 




o 


o 


G 












































3 


X 






ao 


uo 


co 


o 


o 


kD 


CM 


m 


<r 




m 


00 


m 


m 


m 


CO 


o 


<3- 




J-I 










vO 


ON 


CT\ 


i—i 


CM 


co 


o 




<r 


i—i 


CD 


CO 


<r 


cD 




O 


<r 


G 


CO 


00 












i — i 


i — 1 


i—i 


i—i 








1 — 1 


i— i 








i— 1 


00 


o 


00 


00 






































i—l 


PS 


i— 1 


I—l 






































• 


G 


• 


• 






































CL 


o 


CL 


CL 


















































































r— » 


CO 


r—i 


l—l 






































00 


O 


00 


00 






































1 < 


J-i 


i i 


1— .1 








































w 










































1— 1 




i—l 


1— 1 
































-~N 






cu 


£ 


cu 


cu 






























< 


0 






ao 


a 


ao 


ao 
































o 






cu 


G 


cu 


cu 
































ac 






•H 


cu 


•H 


•H 








































PQ 


& 






o 










































o 


o 


o 


o 


o 


o 


O 


o 


o 


m 


m 


m 


m 


m 


m 


333 










CD 


ON 






a\ 


ON 


a>* 






<3- 


<r 




co 


CO 


co 


© 


o 


• 


• 


























i — i 


i — 1 


i—i 


a 


O 


O 


i— 1 


CM 


co 





- 61 - 



Its output is still no better than its input data. Study of 
complete wave spectra are still needed for thorough analysis of an area« 
The program should however make this analysis easier than it has been 
in the paste 

There are still differences between theoretical and experimental 
results as shown by Putnam and Arthur [20], Blue and Johnson [26], and 
Carr and Stelzride [25] , which were mentioned in Section III B, Although 
the agreement between theory and experiment is satisfactory, theory 
does overestimate wave heights by up to 10 percent in the zone beyond a 
breakwater tip. This is probably due to a significant loss in wave 
energy during the process of reflections, (Terry [27]). All calculations 
are based on ideal breakwaters with 100 percent reflection. In practice 
incident waves are only partly reflected and are partly destroyed by 
turbulence (Penney and Price [3]) 0 

A basic assumption is that the wave height is small compared to its 
length and that the wave profile is nearly sinusoidal. Although the ex- 
periments mentioned above indicated that the results are not seriously 
altered by waves of moderate height, design for solitary waves likely 
to be generated during severe storm conditions could be susceptible to 
error. 

Gap solutions have only been experimentally verified by Blue and 
Johnson [26] with a wave angle of approach less than about 60 degrees; 
however, they are the only solutions available so they are actually used 

in practice over the entire 90 degree range. 

For short breakwaters of the ribbon type, with diffraction about 
both ends, the program solutions would not be valid in a region of dif- 
fraction effects from both ends. In practice however, the breakwaters 
are usually of sufficient length to achieve good results by diffracting 
each end independently. 
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IV. PROGRAM SUMMARY AND ENGINEERING FEATURES 



The program developed was written so as to facilitate its use by 
engineers in the analysis of coastal areas. The following summaries 
give highlights of their operation. Details of the main program and 
all subroutines are explained in the appendix. 

(1) Main Program 

The purpose of the main program is to read data and commands, set 
variables, and call on the various subroutines to take appropriate action 
when necessary c The main program is written so as to read problem- 
oriented language input. In this way, the operator uses the same langu- 
age in describing input data and tasks to be performed as he uses in 
normal every day work. When appropriate, a written output of data is 
provided to assure a complete record of the calculations performed. 

An attempt was made to provide most of the commands that an engineer 
would have occasion to use in an analysis of a coastal area. Although 
this program was, of necessity, written for batch processing, it is well 
suited to console input and random interrogation by an operator follow- 
ing the progress of a solution. 

Since much of the refraction/diffraction analysis process is still 
left to engineering judgment, the flexibility was provided for the en- 
gineer to use his judgment whenever it is deemed necessary. A minimum 
of information is set by the program, it depends on an engineer for 
important variables , 
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(2) Refraction Subroutine 



This subroutine traces wave rays across an array of water depths 
until the ray goes onto a beach, goes off an edge of the array, or 
enters a breakwater zone. During this trace, all information concern- 
ing the wave ray which could be of use to an engineer is output, in- 
cluding location, water depth, speed, direction of travel, shoaling 
coefficient, and coefficient of refraction, 

Due to the fact that in many cases refraction analysis is made in 
breakwater areas, this program stops any wave traces upon crossing a 
breakwater or upon entering a breakwater diffraction zone. At these 
points of intersection, all necessary wave ray information is determined 
and output for use by the operator in further analysis « 

(3) Diffraction Subroutine 

This subroutine calculates the coefficient of diffraction and 
direction of wave front travel in a diffraction zone. This may be a 
diffraction zone at the end of a semi- inf inite breakwater, or a diffrac- 
tion zone at a gap between two breakwaters. Given the location of one 
or two breakwaters and necessary data on the incident wave, this progran 
analyzes the diffraction pattern at any specified point and outputs the 
coefficient of diffraction and the wave front direction of travel. 

(4) Shoreline Tracing 

This subroutine will find and trace shorelines across the array of 
water depths. Although this feature is most useful if the program is 
converted to some sort of a graphical output device, it can also be 
helpful when working in an area where tidal fluctuations can have a 
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noticeable effect on the shoreline location, Once a depth array is 
input to the program, presumably based on a "mean-sea level 11 , tidal 
changes can be simulated through the input of water level changes; if 
these changes affect the shoreline, this program will determine its new 
location, 

(5) General Features 

To facilitate maximum adoption of this system by any interested 
users, it was written in the FORTRAN IV programming language. Although 
written for the IBM 360 computer, it can readily be adapted to most other 
hardware simply by changing the alphameric character constants in MAIN 
and by adopting the read/write statements to available devices. 

The program is now written for any water depth array size up to 50 
by 40 grid units. This can be changed up or down as desired, simply by 
changing the DIMENSION cards. No other changes are needed. 

The entire system at present requires about 10,000 words of memory 
on an IBM 360 - Model 40 computer. However for smaller cores, it is 
possible to use the program in portions. For example, during refraction 
calculations, subroutines, RSSHOR which locates the shoreline, and 
RSDIFF which performs diffraction calculations, are not needed and mini- 
mum sized dummy subroutines may be inserted for them. During diffrac- 
tion calculations, only RSDIFF and MAIN are needed. For tracing a shore- 
line, only MAIN and RSSHOR are needed. Although this method of operation 
would cut down on some of the operator flexibility, it can facilitate 
program use on hardware too small to accept the entire package and it 
could enable an operator to trade off unnecessary program for an in- 
creased array size. 
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Vc POSSIBLE SYSTEM EXPANSION AND IMPROVEMENTS 



The programs developed show some of the possible applications of a 
computer by an engineer in the study of .coastal areas. While writing 
these programs it became apparent that many additions and improvements 
could be made to them which would still further facilitate engineering 
analysis* The following are some of these possible expansions. These 
explanations are not intended to be complete analyses but are only to 

i 

point out possible areas of further study. 

(1) Wave Front Tracing 

In many studies, plots of wave fronts can be very useful in 
graphically depicting a situation. The present program, by tracing wave 
rays, will not give any wave front data other than the fact that the 
wave fronts are orthogonal to the ray and simulated fronts could possi- 
bly be added to the rays by an experienced draftsman or engineer. It 
would be possible however to expand the present program to output wave 
front coordinates. The present ray trace program is based on increments 
of equal distance between points. This could be changed to increments 
of equal time. As the iteration process, which is needed to locate 
succeeding points, take place, an average speed would be calculated and 
the interval between points could be based on this average wave speed 
and the given periods of equal time. If all wave rays were started 
along a wave crest, (this ability to start multiple rays could also be 
added to the program) each successive coordinate point output by each 
ray would be a successive crest position. This feature would also have 
the beneficial effect of reducing the increment interval as the ray moves 
into shallower depths. 
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(2) Graphical output 



The refraction/diffraction problem is ideally suited for a 
graphical output. The combination of graphical output and console input 
by an engineer would be an extremely powerful tool in the analysis of 
any coastal area. Early versions of the program were run on an IBM 
1620 connected to a Gerber Plotter. Unfortunately, the program was 
not perfected prior to removal of this facility from M.I.T. 

Although not used by the author, perhaps the most interesting and 
promising output would be accomplished on an oscilloscope. Through a 
combination of console and light pen input to the program, an engineer 
could rapidly analyse diffraction and refraction problems. With a mini- 
mum of effort he could investigate far more possible solutions than 
would be practical by any other method. 

(3) Breakers and Surf 

At present, most predictions of breaker and surf locations along 
a shoaling beach are done empirically . These locations can be very 
important in the analysis of a beach due to the fact that refraction 
does not take place beyond the breaker zone and direction of travel of 
the breaker is very important in determining littoral currents . If 
breaking or surf occurs on a reef at some distance from a shoreline, 
especially at different tidal elevations, this can be very important. 

For some studies therefore, it may be desirable to add to the program a 
means of predicting the breaker location. This could be done by 
determining a desired breaking criteria and adding it as another check 
item to the refraction calculations 0 



- 67 - 



(4) Refraction by Currents 



In some areas, especially at estuaries, one may have to allow for 
refraction caused by currents „ In the past, using graphical methods, 
it was impractical to account for local currents in refraction studies 
and they were either neglected (Dunham [13]) or accounted for by engin- 
eering judgement* A theory for the refraction effects of currents is 
well established and has been published in some detail by Johnson [10] 0 
Due to the possible occurrence of currents in coastal areas resulting 
from long shore currents, currents caused by topography, flow from re- 
gions of high refraction to areas of low refraction, eddy currents, and 
rip currents, it may be advantageous to have an ability to account for 
their effects. This would be a desirable addition to the present 
program. 

(5) Output Summaries 

A suitable means of summarizing output is needed. Various methods 
of portraying refraction for a location, showing the effects of period, 
and direction of wave travel have been used in the past by Johnson [10] , 
and the Beach Erosion Board [9]. As the volume of data increases how- 
ever, most of these methods become confusing. It would appear that this 
is an area requiring considerable common sense on the part of the engin- 
eer to assure that a detailed study with noteworthy results is not lost 
or obscured by a poor presentation. It would be mdst desirable, of 
course, to have a suitable data presentation come directly from the 
computer. 
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(6) Approximations for Diffraction Process 



As noted earlier, diffraction theory, when applied to water waves, 
generally gives conservative results. It was also found that in many 
instances, approximations of the complete theoretical solution brought 
results into closer agreement with the tests than did the complete 
solutions, (Blue and Johnson [26]). The Beach Erosion Board [9] points 
out that becausfe of non-uniform validity of various orders of wave 
theories, the best approximation is not necessarily the highest order, 
LaCombe [30] proposed approximate solutions to the diffraction problem 
which may offer some avenues of approach. It may therefore be desirable 
to adjust the present program to allow for various approximations. 

These approximations can be compared with the complete solution and can 
be checked against experimental data for an empirical comparison. 

(7) Command structure 

To the maximum extent possible, the program was designed to 
facilitate expansion and improvement. Additional commands may be added 
to MAIN simply by adding their appropriate alphameric characters to the 
listing and adding a block of appropriate instructions to the program. 
The existing commands could be improved through the addition of allow- 
able abbreviations and shorter names. Input is well suited to free 
format and this would be a most desirable feature, especially if console 
input is to be used. 

(8) Multiple Waves 

In many instances, different period waves and waves travelling 
from different directions will arrive at a given point at the same 
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time. There are methods for determining the resultant of multiple 
waves as shown by Eagleson, et al [7]. This is not a simple addition 
problem and, in view of its complexity and common occurrence, it would 
be a desirable feature which a design engineer could use. 

(9) Combined Refraction and Diffraction 

All diffraction theory is based on waves of constant speed. In 
practice, breakwaters are usually built in areas where shoaling is an 
important part of the problem and the wave speed variation is important,, 
To date, no combined theory of refraction and diffraction has been 
developed. According to Bretschneider [12], and many other authors, 
the only satisfactory means of studying combined refraction/diffraction 
problems is to refract a wave to a breakwater zone, diffract the wave 
using an average wave speed over an arbitrary number of wave lengths , 
then refract the wave again to the shoreline. Products of refraction 
coefficients and diffraction coefficients are used to determine final 
wave heights. 

Some method of calculating combined refraction/diffraction prob- 
lems would remove these arbitrary evaluations required at present. 



- 70 - 



VI. CONCLUSIONS 



The basic purpose of this study was to develop a computer system 
for use by engineers in the analysis of wave behavior in coastal areas. 

The refraction program developed was based on 3-dimensional con- 
tinuous parabolic interpolation of a depth array. This procedure gives 
continuous curves through any series of data points. Adapted to 3- 
dimensions and used for refraction analysis, it appears to offer 
excellent results. 

The diffraction program is based on a form of the basic Penney- 
Price water wave diffraction methods. The system developed calculates 
diffraction coefficients and directions of wave travel which compare 
closely with manual/graphical plots. 

The entire program package offers a means of analyzing a coastal 
area in whatever depth the engineer desires in but a small fraction of 
the time required by older methods. Hopefully, this will facilitate 
more thorough design than has been economically feasible in the past. 

More testing of the programs is required; hopefully against situa- 
tions where natural conditions are known and where output could be 
compared to actual phenomena. 

Refraction/diffraction analysis is ideally suited for graphical 
plots of the output. With the present program system as a basis, output 
should be adapted to some sort of plotter or oscilloscope where solu- 
tions in progress could be watched and evaluated by an engineer. With 
console input, an engineer watching plots of his solutions would be able 
to try various design schemes for a facility and thoroughly evaluate 
their effects in but a few minutes. 
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Although improvements and further testing are needed, the system 
developed appears to offer considerable promose for better analysis of 
coastal refraction/diffraction problems than has been possible in the 
past ^ 
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VIII APPENDICES 



A. Fresnel Integral Calculations 

B. Program Structure 

Each program includes: a c Listing of Variables 

b. Summary of Operations 

Co Macro-Flow Chart 
d. Program Listing 

1 0 MAIN - Program to read in all commands and data and take or 
call for appropriate action* 

2 C RSRFTN - Subroutine to trace a wave ray through shoaling water < 

3 0 RS1NTP - Subroutine that calculates depth and bottom 

derivatives at any point* 

4o RSBETA - Subroutine which calculates ray separation factor, 
coefficient of refraction and shoaling factor. 

5 0 RSCAFK - Subroutine which calculates wave velocity and 
curvature at any point* 

6. RSDIFF - Subroutine which calculates coefficient of diffrac- 
tion and wave direction of travel in a diffraction 
zone * 

7* RSSHOR - Subroutine which locates and traces shoreline. 

C. Format for Command /Data Input 

D. Refraction Test on Theoretical Beach 

E. Sample Run 
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APPENDIX A 



EVALUATION OF FRESNEL INTEGRALS 



A Fresnel Integral is defined by: 
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using 



Similarly : 
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The following formulae as given by Sparrow [23] can be used to 
approximate C and S. For large values of t, 
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APPENDIX B 



PROGRAM STRUCTURE 
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1. MAIN PROGRAM 



Main program to read data cards and take or call for appropriate 
action. 

PROGRAM NAME: Referred to as MAIN in text. 



Variables stored in Common 



AD. 



1) Angle in radians from -fX direction to wave ray. Angle 
is measured in positive (counterclockwise) direction. 

2) In diffraction, the angle of approach of the incident 
wave ray at the breakwater. 

A in degrees; used for input/output. 



BN, Bl, B2 . . Beta, wave ray separation factor, at point preceeding 
(X,Y) , at (X,Y) and at point succeeding (X,Y) , 
respectively e 

BXAS , BYAS, 

BXBS, BYBS. « . Coordinates of ends of first breakwater, (A-B) on DMAT 
grid. 



BXCS, BYCS , 
BXDS , BYDS o . 



Coordinates of second breakwater (C-D) . If GAP is 
specified it must lie between B and C ends of breakwater. 



C ..... . , Wave velocity in ft/sec. 



CN 



CPPK 



G*T/(2*PI); deep water wave speed. 



C *T/(4*PI); referred to as CK" or K"C in text 
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DCDX, DCDY* . 


. 9C/9X, 9C/9Y in text; partial derivatives of wave speed 

to X and Y, respectively. 


DDDC .... 


. dd/dc in text; derivative of depth to wave speed. 


DMAT .... 


. Matrix of water depths in feet. Positive indicates water 
depth; negative is above shoreline. 


DXY 


. Depth at (X,Y) . ("d" in text) . 


DIX, DIY . . 


o 9d/9X, and 9d/9Y in text; slope of bottom at (X,Y) . 


D2X,D2Y,D2XY 


. 9^d/9X^, 9^d/9Y^, 9^d/9X9Y in text; partial second 

derivatives of bottom depth at (X,Y) . 


FK 


. Wave ray curvature at (X,Y) . Equals inverse of radius 
of curvature. 


G 


e Acceleration due to gravity; 32.2 ft/sec. 


GRID . . . . 


. Size of each DMAT grid unit in feet. 


IXN.IYN . . . 


. Maximum matrix dimensions being used in text. These can 
be less than or equal to dimension statement and are 
used in lieu of changing DIMENSION statement. Can use 
any DIMENSION statement which will fit in memory with 
rest of programs e 


NBW 


o Index for counting breakwaters; NBW=1, if none; 2 if 
one breakwater; 3 if two breakwaters (no gap); 4 if two 
breakwaters and gap between them. 


PI 


o 3.1415927 


PK ..... 


. T/ (4*PI) ; referred to as K' in text. 


PPK ..... 


, (2*PI) / (G*T) ; referred to as K" in text. 



- 81 - 



s. . 



. . Incremental distance, in grid units, between successive 



calculation points along a wave ray. 



T - c Wave period in seconds. 

WLD ...... Length of diffraction zone beyond ends of breakwaters; 

either in feet or number of wave lengths, depending on 
input, 

X,Y o Coordinates of position within DMAT which are being 

studied . 



Variables used by MAIN 



DMMY ....oo 
DX , DY , DM . . . . 



I , J • . • . • • 

IX, IY 

1X1 ,1X2 , 

IY1 , IY2 . . . . 
IWD (— ) . . . . 

IWD1 ,IWD2 . . . . 
IXP1 , IYP1 . . o o 



Dummy variable, no significance. 

Used in establishing analytical depth matrices; DM is 
the depth at DMAT (1X1, IY1) ; DX is the increase in 
depth per X grid unit; DY is the increase in depth per 
Y grid unit. 

"Do Loop" indexes. 

Grid coordinates on DMAT (IX, IY) array. 



Limits of operations being performed on DMAT. 

Listing of decimal equivalents of "A" format characters 
for IBM 360. 

"A" format characters read off command portion of input. 
Convenient totals IX + 1; IY + 1. 
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LL,LL2 



Index used to check through IWD(LL) array. 



LLS . . . . o . Value of LL read on last command. Needed if present 

command is a blank which indicates repeat of last 
command. 

NN Index used to indicate program status when calling 

RSCAFK: NN=1, cako depth and only C; NN=2 , depth is 
known, calc. C and FK. 

NWL. o .... o Index to indicate means of measuring length of 

diffraction zone. (See WLD) . 

NWL=1 if measured in wave lengths. 

NWL=2 if measured in feet* 

N1,N2,F1,F4 . . Fixed and floating point data values read off command 

card. These will be set to actual variables after command 
is decoded. 



Summary of Program 

The purpose of MAIN is to read all command and data cards, then, 
either within itself or by calling subroutines, take whatever action is 
required. Its principal functions are bookkeeping, establishing variables 
for use by subroutines, and calling subroutines to take whatever action 
has been called for by a command. In some instances, variables are pre- 
determined; in some, a value is assumed unless changed by the operator; 
and in others, operator input is essential. The proper format of all 
commands is shown in Appendix C. 



- 83 - 



Command /Data Input 



See Appendix C for proper formate Underscore indicates Alphameric 
characters which are used in decoding command and which must occur in 
columns 1 and 10* Abbreviations which do not change underscored charac- 
ters are allowed. 

Required input for Refraction Calculations : 

Either an analytical or topographic array or a combination of the two 
must be established. 



ANALYTICAL MATRIX 


1X1 


1X2 


DX 


DM 




IYl 


IY2 


DY 





Two cards are required, A plane surface is established which extends 
from 1X1 to 1X2 with an increase in depth of DX per X grid unit; and ex- 
tends from IY1 to 1Y2 with increases of DY per Y grid unit. DM is the 
depth of (1X1, IY1) . All depths are in feet. 

LIMITS 0F_ TOPO 1X1 1X2 

IYl IY2 

Followed by: 

TOPO IN FEET FOLLOWS 
(or) 

TOPO IN FATH FOLLOWS 
Followed by: 

Data in 10 F6.4 Format 
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This is a series of commands used in reading natural beaches into 
array which extend from 1X1 to 1X2 and IY1 to IY2. If topography is in 
fathoms, the program immediately converts it to feeto 

Data is read one X row at a time, 10 items per card« 

WATER LEVEL FI 

At any time in the program, the water level may be changed to 
simulate tidal fluctuations . This card increases all water depths by Flo 
Lower water levels can be specified by negative Flo 

RAY DATA JL NTC T AD X Y 

This card establishes the starting point of a wave ray at (X,Y) with 
an angle AD degrees from the + X axis c NTC is any integer control number 
desired by the operator for ray identification. T is the wave period in 
seconds . 

RAY DATA 2 B1 BN 

This is an optional card which is used only if a ray separation 
factor, other than 1.00, is desired at the start of a wave ray. Bl is 
the Beta value at (X,Y) , BN is the Beta value at an imaginary point 
(XN,YN) preceeding (X,Y) by a distance S along the ray. 

INCREMENT S 

This is an optional card which can change the initial increment 
used in tracing the wave ray from 4.0 to 2*0, 1.0 or 0,5 grid units. 

If not specified, a value of 4.0 is used. Any interval specified is a 
maximum and the program will automatically change to lesser values as 
the ray curvature increases «> 
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TRACE RAY_ 

This traces wave ray specified by RAY DATA 1 and RAY DATA 2 until it 
goes off edge of array, runs onto beach, or strikes breakwater zone. 
Program will output successive X, Y coordinates of ray; wave speed, C; 
ray curvature, FK; ray separation factor, BETA; coefficient of refrac- 
tion, FKD; The Shoaling Coefficient, D; Ray Angle, AD; and the current 
increment used in tracing ray, 3„ 

Commands required if breakwaters are desired on array : 

Lack of input assumes no breakwaters; one or two breakwaters may 
be specified. 

BREAKWATER BXAS BYAS BXBS BYBS 

Establishes a breakwater (A-B) from (BXAS, BYAS) to (BXBS, BYBS). 

BREAKWATER BXCS BYCS BXDS BYDS 

Establishes a second breakwater (C-D) from (BXCS, BYCS) to 
(BXDS, BYDS). The first breakwater specified is automatically designated 
A-B; second is C-D* 

GAP 

Specifies a breakwater gap from (BXBS, BYBS) to (BXCS, BYCS). If 
not specified, the two breakwaters are considered independent. This card 
can only be used if two breakwaters have previously been input to program* 

SIZE GRID_ IXN IYN GRID 

Sets limits of array with X from 1 to IXN and Y from 1 to IYN. If 
not specified, 50 by 40 is assumed. GRID is the size in feet of each 
grid unit of array. 
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Optional Breakwater Commands : 



WAVE LENGTHS DIFFR WLD 

This specifies the effective distance, in wave lengths, of a break- 
water beyond its tip into the diffraction zone. The wave ray will be 
traced until it either strikes a breakwater or passes within this number 
of wave lengths of a breakwater tip as measured along an orthogonal to the 
ray. Upon striking a breakwater or breakwater zone, ray values at the 
point of intersection are calculated« 

LENGTH DIFFRACTED WLD 

Similar to WAVE LENGTHS DIFFR except that WLD, in this case, is 
expressed in feet. If neither diffraction distance is specified, a value 
of 5.0 wave lengths is assumed. 

ELIMINATE BWS 

This command removes all breakwaters and gaps which have been specified. 
New breakwaters may afterwards be specified if desired. 

Input required for diffraction calculations : 

The location of one breakwater (A-B) , or two breakwaters (A-B) and 
(C-D) must be specified as called for under refraction commands. 

DIFFRACT A 
DIFFRACT B_ 

DIFFRACT C 
DIFFRACT D 
DIFFRACT GAP 

Any one of these commands specifies that diffraction calculations 
should take place about the designated breakwater tip. One breakwater 
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must have already been input if A or B is specified; two if C or D is 
specif ied. 

If GAP is specified, an imaginary gap normal to the wave ray is 
established between points B and C. B and C must already be in the pro- 
gram 0 GAP also requires that a wave ray angle, A, be in the program. 

If one remains from, refraction calculations, it will be used; if not, or 
if a change is desired, a RAY DATA 1 card must be input prior to DIFFRACT 
GAP. 

Any changes in A must be followed by a new DIFFRACT GAP Command. 
DIFFRo COORD X Y 

Calculates the coefficient of diffraction and direction of wave 
travel at (X, Y) . In addition to having already specified which break- 
water tip is to be diffracted through a DIFFRACT command; a wave speed, 

C; wave ray angle, A; period, T; and grid size, GRID; must already have 
been specified e GRID must be input through a SIZE GRID Command. If 
present, T and A may be left unchanged from prior operations or may be 
specified by a RAY DATA 1 card* If present, C may be left unchanged from 
a prior portion of the program or must be specified by a CALC WAVE SPEED 
or by a WAVE SPEED Command. 

Other Commands : 

WAVE SPEED C 

Sets C equal to a specified value. 

CALC WAVE__SPEED X Y 

Calculates and sets C equal to the wave speed at (X, Y) . 
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LOCATE SHORELINE 

Based on depth data, which must already have been specified, this 
command locates and traces the shoreline across an array* Shoreline 
coordinates are output*, 

_(blank) _ N1 N2 FI F2 F3 F4 

Indicates repeat of preceeding command with new data values* 

ALL DONE __ 

Calls for end of the program* 
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MACRO FLOW CHART OF MAIN 



Bookkeeping and 
setting of constants 




n n n 



CALL FOR APPROPRIATE ACTION. 
CAMBRIDGE* MASS* 1967 



C MAIN PROGRAM TO READ IN DATA AND 

C PROGRAMMED BY R D SMART AT M I T 

DIMENSION DMAT ( 50 *40 ) , IWD( 13 ) 

COMMON DMAT * I XN» I YN »P I »G»GR ID*X *Y *DXY *D1X *D1Y »D2X »D2Y, D2XY* 

1 A*AD*T,FK*BN*B1 , 32 * C * PK » PPK , CPP K , CN » DDDC » DCDX » DCDY , 

2 S* NBW*BXAS *BYAS*BXBS*BVpS»BXCS*BYCS*BXDS*BYDS*WLD 
199 FORMAT ( A 1 ♦ 8X » A1 * 1 OX , 2 I 5 ♦ 5F 1 0 . 4 ) 

( 10F6.4 ) 

(19H TILT INPUT ERROR * 4HWD 1 = * A 1 » 6H WD2=*A1) 

(20H DEPTH ADJUSTMENT IS*F5.3* 3HFT.) 

(10H ALL DONF. ) 

(18H WAVF SPEED AT X= *F5.2* 3H Y= *F5.2* 4H IS * 
16H FT/S* DEPTH IS ,F6.3*4H FT. ) 

WAVE SPEED IS SET AT ,F6.2*5H FT/S ) 



FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
1 F6.2, 

193 FORMAT (22H 



198 

197 

196 

195 

194 



» I 2 * 1 1H PERIOD = 
4HDEG . , 2HX = ♦ F5.2* 



» F5 • 1 » 4HSEC » * 
2 H Y = » F5.2) 



,F6.2*4H TO * F6 . 2 , 



192 FORMAT (12H TEST CASE 
1 8H ANGL E= * F5.1* 

191 FORMAT (20X»2I5*F10.5) 

(17H BREAKWATER FROM ,F6.2*2H. 

2H, *F6.2 ) 

( 5H GAP ) 

( 24H ELIMINATED BREAKWATERS ) 

( 26H INITIAL RAY INCREMENT IS , F 5 - 2 ) 

(15H DIFFR DIST IS ,F8.2, 6H FEET ) 

(15H DIFFR DIST IS *F6.2*14H WAVE LENGTHS ) 
(25H ANALYTICAL DATA FROM X= »I4*4H TO ,14, 
8H AND Y = ,14 ,4H TO ,14 ) 



190 FORMAT 

1 

189 FORMAT 
FORMAT 
FORMAT 

format 

FORMAT 
FORMAT 



188 

187 

186 

185 

184 



1 

183 FORMAT 
1 

182 FORMAT 
1 



(9H GR ID IS * 14, 4H BY 
15H FOOT SQUARES. ) 
(19H TOPO DATA FROM X= 
14, 4H TO ,14) 



, 14, 6H WITH 



r 8 • 2 



, 14, 4H TO , 14, 8H AND Y = 



M I SCFL L AN EOUS BOOKKEEPING 



I WD ( 1 ) = -1052753856 
I WD ( 2 ) = -1035976640 

I WD ( 3 ) = -10191 99424 

I WD ( 4 ) = -1002422208 
I WD ( 5 ) = -952090560 

I WD ( 6 ) = 1077952576 

I WD ( 7 ) = -650100672 
I WD ( 8 ) = -918536128 

I WD ( 9 ) = -482328512 
I WD ( 10)= -431996864 
I WD (11) = -750763968 
I WD( 12 ) = -499105728 
I WD (13) = -985644992 
PI = 3.1415927 
G = 32.2 

C ASSUMED CONSTANTS UNTIL RESFT BY OPERATOR 

DO 50 IX = 1,50 
DO 50 IY = 1,40 
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50 DM AT ( IX , I Y ) = 

N8W = 1 
B 1 = 1 • 

R2 = 1 . 

S = 4. 

NWL = 1 
WLD = 5. 

IXN = 50 
IYN = 40 

READ IN ONE COMMAND/DATA CARD 

100 READ (5,199) I WD1 » I WD2 > N 1 > N2 » F 1 , F2 , F 3 » F4 
LLS = LL 

CHECK IWD1 AGAINST ALLOWABLE INPUT CHAR ACTERS . . ( A FORMAT) 

DO 104 LL = 1,13 
IF (IWD(LL) - IWD1) 104,106,104 
104 CONTINUE 
108 WRITE (7,197) 

GO TO 138 

IF NO COMMAND IS GIVEN (BLANK), ASSUME REPEAT OF OLD COMMAND 
150 LL = LLS 

GO TO LOCATION SPECIFIED BY COL. 1 

106 GO TO ( 138, 1?0 , 137 , 110, 126, 1 50, 128 , 134, 114, 116, 140, 1 12, 142 ) ,LL 

READ A DIFFRACTION COMMAND 

1 10 DO 1111 LL2 = 1,7 

IF ( I WD2- IWD ( LL2 ) ) 1111,1112,1111 

1111 CONTINUE 
WRITE (7,197) 

GO TO 138 

1112 CALL RSDIFF ( LL2 , I WD2 »F1 , F2 ) 

GO TO 100 

READ IN GRID SIZE 

112 IXN = Ml 
IYN = N 2 
GRID = FI 

WRITE (7,183) IXN, IYN, GRID 
GO TO 100 

114 IF ( I WD2 - 1077952576) 1141,136,1141 
READ IN GRID DATA 
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1141 READ ( 5 ,198 ) ( ( DMAT ( I , J) , I = 1X1 , 1X2 ) » J= IY1 » I Y2 ) 
C WHEN NECESSARY CONVERT FROM FATHOMS TO FEET 
I F ( I WD2 + 1052753856) 100,115,100 

1 15 DO 1151 1 = 1X1 ,1X2 

DO 1151 J= I Y 1 , I Y2 

1151 DMAT (I,J) = 6 # *DM A T ( I , J ) 

GO TO 100 

116 IF ( IWD2 + 1002422208) 1161,1371,1161 

1161 I F ( I WD2 + 482328512) 1162,1163,1162 

SET DIFFRACTION IN WAVF LFNGTHS 

1163 WLD = FI 
NWL = 1 

WRITE (7,185) FI 
GO TO 100 

ADJUST WATER LEVEL 

1162 WRITE (7,196) FI 
DO 118 1 = 1, IXN 
DO 118 J=1 , IYN 

118 DMAT ( I , J ) = DMAT(I»J) + FI 
GO TO 100 

READ IN BREAKWATER COORDINATES (2 MAX.) 

120 GO TO ( 122,124,108) , NBW 
FIRST BREAKWATER 
122 BX AS = FI 
BY AS= F 2 
BX BS= F 3 
BYBS=F4 
NBW = 2 

WRITE (7,190) FI ,F2,F3,F4 
GO TO 100 
SECOND BREAKWATER 
124 BXCS = FI 
BY CS = F 2 
BXDS = F 3 
BYDS = F 4 
NBW = 3 

WRITE (7,190) F1»F2»F3»F4 
GO TO 100 

SPECIFY GAP RFTWPFN 2 BRFAKWATFRS 

126 I F ( NBW -3) 108,127,108 

127 NRW = 4 
WRITE (7,189) 
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GO TO 100 
RFAO IN RAY DATA 

128 IFUWD2 + 247447488 ) 1 30* 132. l’O 
132 NTC = N 1 
T = FI 
AO = F 2 

A = AD*P 1/180. 

X = F 3 

Y = F4 

GO TO 100 
130 B 1 = FI 
B2 = F 2 
GO TO 100 

CALC ANALYTICAL MATRIX 

138 IF ( I WD2 + 750763968) 1381,1382,1381 

1382 1X1 = N 1 
1X2 = N 2 
DX = FI 
DM = F 2 

READ (5,191) I Y 1 , I Y 2 , DY 
DMAT ( I X 1 , I Y 1 ) = DM 

IXP1 =1X1+1 
IYP1 = IY1 + 1 

DO 1385 IX = IXP1 , 1X2 

1385 DMAT ( I X , I Y 1 ) = DMAT ( I X - 1 , I Y 1 ) + DX 

DO 1387 IY = I YP1 , I Y2 

DMAT ( I X 1 , I Y ) = DMAT( I X 1 , IY - ] ) + DY 

DO 1387 IX = I XP1 , 1X2 
1 387 DMAT (IX, IY) = DMAT (IX - 1 , I Y ) + DX 

WRITE (7,184) 1X1 » 1X2 » I Y1 , IY2 

GO TO 100 

ALL DONE QUIT 

1381 WRITE (7,195) 

CALL EXIT 

CALCULATF WAVE SPFFD 

137 NN = 1 
X = FI 

Y = F2 

CALL RSCAFK. (NN) 

WRITE (7,194) X , Y , C * DXY 
GO TO 100 

PRESET WAVE SPFFD (FOR DIFFRACTION CALCULATIONS) 
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c 

1371 C = FI 

WRITE (7*19?) C 
GO TO 100 

TRACE RAY 

136 WRITE (7*192) NTC*T*AD*X*Y 
CALL RSRFTN (NWL) 

RESET NORMAL INITIAL CONDITIONS 
B 1 = 1. 

B2 = 1. 

S = 4. 

GO TO 100 

CHANGE NORMAL TRACE INTERVAL FOR RAY CALCULATIONS. 

134 S = FI 

WRITE (7*187) FI 
GO TO 100 

140 IF ( IWD2-1077952576 ) 1401,1402*1401 

1401 IF ( I WD2 + 968867776) 1403*1404*1403 

SET DIFFRACTION DISTANCE IN FEET 

1404 NWL = 2 
WLD = FI 

WRITE (7,186) FI 
GO TO 100 

TRACE SHORELINE 

1403 CALL RSSHOR ( DMMY ) 

GO TO 100 

READ IN LIMITS OF TOPO DATA 

1402 1X1 = N 1 
1X2 = N 2 

READ (5*191) I Y 1 , I Y2 

WRITE (7*182) 1X1 , 1X2 * IY1 , IY2 

GO TO 100 

ELIMINATE BREAKWATERS 

142 NBW = 1 

WRITE (7,188) 

GO TO 100 
END 
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2o SUBROUTINE RSRFTN 



Subroutine which traces a wave ray through shoaling water . 

SUBROUTINE NAME: RSRFTN 

Variables Used in Program ; 

COMMON Variables are listed under MAIN* 

* 

ABAR . . . . « Ray Angle from (XN, YN) to (X, Y) ; an average of the 
angles at the two points c 

AN, AMI .... Last 2 preceeding values of A. Used to interpolate 
intermediate value between AN and A. 

BBW1 Y intercept of Breakwater No. A-B. 

BBW2 Y intercept of Breakwater No. C-D. 

BBW3 Y intercept of breakwater gap. 

BL Y intercept of breakwater component. 

BR Y intercept of the ray. 

D Shoaling Coefficient; D in text 

DELS In case of intersection, the distance from (XN , YL) as 

compared to the original S distance. 

DELTA Change in Ray Angle, A, from (XN, YN) to (X, Y) . 

DMMY Dummy; no significance. 

DN 2*(Deep water wave length) ; used to estimate deepest 

reasonable depth at which waves will feel bottom. At 
depths greater than DN, it is assumed the wave ray 
curvature is zero. 
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DN2 



0*5 (deep water wave length) ; depth at which waves are 
substantially affected by bottom. Between DN and DN2 , 
only one wave ray curvature calculation is made. At 
depths less than DN2 , up to 10 iterations are used to 
arrive at a final curvature, 

FKBAR Value of FK used to estimate curvature from (XN, YN) to 

(X, Y) ; usually it is an "average” value. 

FKD ...... Coefficient of Refraction; K, text. 

d 

FKL ...... The last value of FK; used to check for convergence at a 

given (X, Y) point. 

FKN, FKNMI . . Last 2 preceeding values of FK. Used to interpolate 
average value between FKN and FK. 



FLEN Length of diffraction zone in grid units. 

FLM Slope of the breakwater component or zone being -checked 

for intersection with the ray. 



FLMR Slope of the ray from (XN, YN) to (X, Y) . 



FLMR2 Inverse slope of FLMR. 

FMBW1 Slope of Breakwater No, A-B. 

FMBW2 Slope of Breakwater No. C-D. 

FMBW3 Slope (x/y) of breakwater gap. 

N Number of intervals moved along ray. Ray stops at 

N=100 as a reasonable maximum. 

NDXY Index used to determine location in program prior to 

breakwater intersection test; 1 if in refraction zone, 
2 if in deep water. 
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NFKI 



Counter for the number of FK iterations at a given point. 
If not satisfactory convergence occurs within 10 itera- 



tions, an average value is used and the ray continues. 



NH Intersection test index 0 

NN Index for RSCAFK to indicate the program status when 

called; NN=1 indicates need to have DXY calculated by 
RSINTP and only C calculated by RSCAFK; NN=2 indicates 



DXY is known, C and FK are desired. 

NP Index for RSBETA: NP=1, at N=l; NP=2 , at N greater 

than 1; NP=3, at a breakwater intersection. 

NT Index to assure that, when necessary, at least 2 FK 

iterations are taken. 



NWL Indication from MAIN as to whether WLD is in feet 

(NWL=2) or WLD is in wave lengths (NWL=1) . 

SOLD Old value of S, 



XDIST , YDIST . X and Y components of FLEN. 

XL, YL . . a . X and Y intersection coordinates of the ray with the 

breakwater component . 

XN, YN . . . c Preceeding (X, Y) point. 

XI, X2, 

Yl, Y2 .... X and Y limits of breakwater component. 
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Functions used in Program 



BNEWF (AX, BX, CX, DX) Given 3 values, AX, BX, CX, function uses 

parabolic interpolation to calculate value of 
a point which is "DX" fraction of way from BX 
towards CX. 

Summary of Program 

RSRFTN is called by MAIN after a DMAT array is established and 
after initial ray values have been assigned, including the starting 
point, (X, Y) , period, T, bearing angle, A, and the initial wave ray 
separation factors, BN and Bl. RSRFTN uses iteration procedures to cal- 
culate the path of a ray through a DMAT array of varying depths. RSINTP 
is called to calculate water depth and bottom derivatives through the 
use of 3-dimensional continuous parabolic interpolation. RSCAFK is 
called to calculate wave speed, C, and curvature, FK, at given points; 
RSBETA is called to calculate the wave ray separation factor, Beta, the 
coefficient of refraction, , and the shoaling coefficient, D. The 
workings of these subroutines are described in succeeding sections. 

An iteration process is used to locate (X, Y) from a preceeding 
(XN, YN) point. Point (XN, YN) , the ray curvature at that point, FKN, 
and the ray angle, AN, are used to estimate a location, (X, Y) , of the 
ray some incremental distance, S, ahead. RSINTP is then called to 
determine bottom depth, DXY , and the bottom derivatives at the new (X, Y) . 
These data are then used by RSCAFK to calculate C and FK at (X, Y) . This 
new FK is then averaged with FKN to give FKBAR. This new average curva- 
ture is then used to make a revised estimate of (X, Y) . At this new 
(X, Y) , the curvature is again calculated, a new average curvature is 
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determined and the (X, Y) estimate is again revised. This estimating is 
continued until FK does not change more than 0.0005 from the FK at the 
preceeding (X, Y) estimate. Once (X, Y) is established, a check is made 
to determine if the ray from (XN, YN) to (X, Y) intersected any break- 
waters or passed within any breakwater diffraction zones. If no inter- 
section occurred, RSBETA is called to calculate Beta, , and D. Upon 
return, the ray values at (X, Y) are printed. 

The ray is then incremented to point (X, Y) which becomes (XN, YN) 
and the process begins again. 

If at any point, DXY is negative, it is assumed a shoreline has been 
crossed. If the ray passes within 2 grid units of the array edges, the 
ray is stopped. 

In water deeper than 2 wave lengths, FK is assumed to be zero and 
Beta remains constant while the ray is incremented along a straight line. 

In depths between 0.5 and 2 wave lengths, only one calculation of 
FK is made and no iterations are used. 

As the ray is traced shoreward, S is reduced from a maximum value of 
4.0 to 2,0 as soon as FK is greater than 0.0001; again from 2.0 to 1.0 
when FK is greater than 0.0005, and finally from 1,0 to 0,5 when FK is 
greater than 0,01. S never exceeds its old value however. Initial S 
values shorter than 4,0 may be specified; this will not change the 
reduction process. 

Slope intercepts are used to determine intersection of the ray with 
any breakwater or breakwater diffraction zone (see MAIN) , Breakwater 
and gap slopes and Y-intercepts are determined at the start of RSRFTN 
for use later as needed. Because diffraction beyond breakwater tips is 
considered to be orthogonal to the ray angle, these slopes and Y-inter- 
cepts must be computed individually. If intersection of a ray and 
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breakwater component occurs, ray values at the point of intersection are 
calculated directly or by using a quadratic approximation and the ray is 
stopped. 

A few safeguards are included in the program to prevent excessive 
running time in case of error . 

N is limited to 100 as a reasonable maximum number of iterations for 
a ray to require before terminating at a shoreline or an edge. 

NFK1 is limited to 10 as the maximum number of iterations allowed 
for convergence of FK values at any one point. Upon reaching 10 iterations, 
a "curvature approximation” message is printed, an average FK value is 
used, and the ray continues. 
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SUBROUTINE RSRFTN ( NWL ) 

C SUBROUTINE WHICH WILL TRACE A WAVE RAY THROUGH SHOALING WATER. 
C BASED ON 3-DIMFNSIONAL PARABOLIC INTERPOLATION. 

DIMENSION DMA T ( 50 » 40 ) 

COMMON DMAT.IXN.IYN.PI » G ♦ GR I D ’ X » Y , DX Y » DlX » D 1 Y »D2X»D2Y»D2XY » 

1 A , AD , T .FK.BN.Bl , B2 , C » P K , PPK , CPPK * CN , DDDC , DCDX , DCD Y * 

2 S.NBW.BXAS tBYAS.BXBS,BYRS.BXCS»BYCS,BXDS.BYDS»WLD 



BOOKKEEPING AND SETTING OF CONSTANTS. 



FUNCTION TO CALC INTERMEDIATE VALUE IN A SERIES 

BNFWF( AX.BX.CX.DX) = BX+ ( CX-AX ) *DX/2 . + ( ( CX-2.*BX+AX ) *DX **2 ) /2 • 
WRITE (7,992) 

992 FORMAT (38X, 7H COFF , 5HSH0AL ) 

WRITE (7.993) 

993 FORMAT (2X.3H N ,7H X , 7H Y 



1 6H BETA ,7H RFT N 
1 PK = T / ( 4 • *P I ) 

PPK = 2.*PI/(G*T) 
CN = G*T/(2.*PI) 

C = CN 

DN = ?.0*CN*T 
DN2 = 0.5*CN*T 
NDXY = ] 

BN = B 1 
AN = A 
AM] = AN 
XN = X 
YN = Y 
N = 1 
NP = 1 
NT = 1 



3H C0EF.8H CURV 



»7H 
, 5H 



»6H C 
,6H DPTH ) 



C 

C SET UP CONSTANTS FOR BREAKWATER INTERSECTION EQUATIONS 
C WATCH OUT FOR INFINITE SLOPES 

C NBW = NO. OF BREAKWATERS... 1 = NONE.... 2 = 1 BW.... 

C 3=2 BW.... 4=2 BW+ GAP 

C 

GO TO ( 10,103.102.101 ), NBW 

101 IF ( ABS< BXCS-RXBS ) - O.OOOOl) ] 00 1 , 1 00 1 , 1 00 2 

1001 BXCS = BXC^ + 0.0001 

1002 FMBW3 = ( BYCS-BYBS ) / ( BXCS-BXBS) 

BBW3 = BYBS - FM0W3*BXBS 

102 IF ( ABS( PXDS-BXCS ) - 0.00001) 1 003 , 1 003 , 1 004 

1003 BXDS = BXDS + 0.0001 

1004 FMBW2 = ( BYDS-BYCS ) / ( BXDS-BXCS) 

BBW2 = BY CS - FMRW2*BXCS 

103 IF ( ABS ( RXRS - RX A S ) - O.OOOOl) 1005,1005,1006 

1005 BX AS = BX AS + 0.0001 

] 0 0 6 FMBW1 = ( BYBS-RYAS ) / ( BXBS-BXAS) 

BBW 1 = BYAS - FMBW1 *BXAS 
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FIND DEPTH AND BOTTOM SLOPES AT X,Y 

10 CALL RSINTP ( DMMY ) 

CHECK FOR SH0R c L I NF 
IF(DXY) 10C, 100, n 

CHFCK TO SEE IF DEFP WATER DEEP = 2*WL 

IS IF (DN-DXY) 12*12*50 

IF DEEP WATER. .NO TRIAL AND ERROR GO ON WITH FK = 0 > DEL BETA = 0 

12 C = CN 
FKN = FK 
FK = 0. 

BN = B 1 
B 1 = B 2 
FKD = 1. 

D*l. 

CHECK FOR BRFAKWATER 

NDXY = 2 
GO TO 90 

19 NDXY = 1 

ITERATIONS HAVF STOPPED AT THIS POINT 
WRITF OUT RESULTS 

20 AD = A * 180.0/ PI 

WRITE (7,998) N * X , Y , AD * C * B 1 * FKD * D * F K , S , DX Y 
998 FORMAT ( 1X,I3,3F7.3»F6.2,F6.3,F7.3»F5.2»F8.4,F5.2,F6.2 ) 

CHFCK FOR EDGFS 

22 IF ( ( X-2 . ) *( FLOAT ( I XN ) - 1. - X)) 100,100,23 

23 IF ( ( Y — 2 • ) * ( F L 0 A T ( IYN) - 1. - Y)) 100,100,24 
C SAFETY STOP FOR RAY AT N = 100 

24 IF C N— 100) 30,30,1 00 
C 

C S CHECK CHANGE S VALUE DEPENDING ON CURVATURE OR RAY. 

C IF FK L.T. 0.0001 USE S = 4.0 GRID UNITS. ..THEN USE S = 2.0 UNTIL 

C FK G. T . 0.005, THEN USE S = 1.0 UNTIL FK G. T. 0.01, 

C THEN USE S = 0.5 
C DO NOT FXCEED OLD S VALUE 

C UPON CHANGING S ADJUST B2 TO INTERMEDIATE VALUE. 

C 

30 SOLD = S 

IE (S-3.999) 202,202,210 

202 IF (S- 1.999) 204,204,220 

204 IF (S-0.999) 300,300,230 
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210 


I F 


( ABS ( FK ) 


-0.0001 ) 


300, 000, 220 




220 


I F 


( ABS(FK) 


- 0.005 ) 


222,222,230 




222 


IF 


(SOLD - 2. 


QOl ) 300 


,300 ,224 




224 


S = 


2. 








225 


B2 


= BN r WF ( BN 


»B1 . R2 .0 


.5 ) 






GO 


TO 300 








2 30 


I F 


( ABS( FK ) 


- 0.01 ) 


232,232,240 




232 


IF 


( ( SOLD - 1 


.001 ) - 


0.0001 ) 300,300 


,234 


234 


S = 


1 .c 










IF 


( (SOLD -2. 


001 ) - 0 


.0001) 225,225, 


236 


236 


B 2 


= BNEWF ( BN 


.81 . B2 .0 


.25 ) 






GO 


TO 300 








240 


S = 


0.5 










IF 


( ( SOLD -1 . 


001 ) - 0 


.0001) 225,225, 


244 


244 


IF 


( ( SOLD - 2 


. 001 ) - 


0.0001) 236,236 


,248 


248 


82 


= BNEWF (BN 


.B1.B2.0 


.125 ) 





c 

C INCREMENT TO NEXT POINT BASED ON BEST FKBAR ESTIMATE. 

C WHEN INITIALLY INCREMENTING. SAVE LAST POINT AND ASSOCIATED 
C VALUES AS A JUMPING OFF POINT FOR NEXT POINT. 

C 

300 AMI = AN 
AN = A 
XN = X 
YN = Y 
FKNM1 = FKN 
FKN = FK 
NT = 1 
N = N+l 
F< BAR = FK 

GO ON TO NEXT X.Y POINT BASED ON BEST FKBAR ESTIMATE. 

AO DELTA = S * FKBAR 
A = AN + DELTA 
ABAR = AN + delta/?. 

X = XN + S*COS(ABAR) 

Y = YN + S*SIN(ABAR) 

GO TO 10 

CALCULATE WAVE SPEED AND CURVATURE AT POINT (X.Y) 

50 IF (NT-1) 52.52.80 

C MAKE FIRST ESTIMATE OF FK 
52 NFK I = 1 
NT = 2 

C CALC WAVE SPEED AND CURVATURE (C AND FK) 

54 NN = 2 

CALL RSCAFK ( NN ) 

C N = 1 POINT ONLY REQUIRES AN FK VALUE. 

IF (N-l ) 920,920.541 
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541 GO TO (70»7 n »82)»NT 
70 FKL = FK 

FKBAR = BNFWF ( F K N M 1 »FKN»FK»0.5) 

MAKE JUST ONE PASS IN )) INTERMEDIATE) ) DEPTHS 
INTERMEDIATE = L • T « 2WL » OR G.T. 0.5WL) 

IF ( DXY - DN2 ) 40,90,90 

IN SHALLOW WATFR RFVISE FK ESTIMATES UNTIL TWO AGREE WITHIN 
0.0005. . .MI NIMUM OF TWO ESTIMATES ARE NEEDED. 



80 NT = 3 

NFKI = NFK I +1 
GO TO 54 

CHECK FOR CONVERGENCE.. ..DONT TRY MORE THAN 10 TIMES 
82 IF ( (ABS(FKL-FK) ) -0.0005) 90,00,84 
84 IF (NFKI - 10) 70,70,86 
86 WRITE ( 7,995 ) 

995 FORMAT (24H CURVATURE APPROXIMATED ) 

FK = (FKL + F K ) / 2 » 

CHECK FOR INTERSECTION WITH BREAKWATER OR BREAKWATER ZONE OF 
OF INFLUENCF. .. . IF THERE ARE ANY BREAKWATERS... 

90 IF ( NBW - 2) 92,255,256 

255 IF (ABS(X-XN) -0.00001) 2551,2551,256 
2551 IF ( ABS ( Y-YN ) -0.000C1 192*92,256 

C CHECK FOR BREAK WATFR NO.A-B 

256 FLMR = ( Y-YN ) / ( X-XN ) 

BR = YN - FLMR*XN 
FLM = FMBW1 

BL = BBW1 
XI = BX AS 
X 2 = BX BS 
NH = 1 
GO TO 292 

C CHECK FOR DIFFRACTION ZONF OFF A END. 

257 IF (GRID-0.0001) 2574,2574,2575 

2574 WR ITE ( 7,299 ) 

299 FORMAT ( 52H PLEASE TELL ME MY GRID SIZE, OR I WONT WORK FOR YOU.) 
GO TO 100 

2575 GO TO ( 25 7 1 , 2 572 ) , NWL 

2571 FLFN = WLD * C * T / GRID 
GO TO 2573 

2572 FLFN = WLD / GRID 

2573 IF ( ABS ( FLMR ) -0.0001 ) 2578,2678,2576 
2578 XDIST = 0. 

GO TO 2577 

2576 FLMR 2 = -l./FLMR 

XDIST = FLEN/( l.+FLMR2**2)**0.5 
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2577 X2 = BX AS - XDIST 

YD 1ST = FLFN/ ( 1 .+FLMR**2 ) **0. 5 
IF ( BYRS- P YAS ) 262*262*260 
260 Y0IST = -YDIST 
262 Y 2 = BYAS + YDIST 
XI = BX AS 
Y 1 = BYAS 
NH = 2 
GO TO 290 

C CHECK TO SEE HOW MANY BREAKWATERS THERE ARE. 

2 6 A GO TO (92*266*266*272 ) ,NBW 

C CHECK FOR DIFFRACTION ZONE OFF END B 
266 X 2 = BX BS + XDIST 
Y 2 = BY BS - YDIST 
XI = BX BS 
Y 1 = BYBS 
NH= 3 
GO TO 290 

C CHECK FOR DIFFRACTION ZONF OFF FND C. 

268 IF (NBW-2 )92 *92*270 
270 X 2 = BXCS - XDIST 

IF ( ( BYDS-BYCS )*( BYBS-BYAS) ) 2701,2702*2702 

2701 YDIST = -YDIST 

2702 Y 2 = BY CS + YDIST 
XI = BXCS 

Y 1 = BYC5 
NH = A 
GO TO 290 

C CHECK FOR INTERSECTION WITH GAP. 

272 FLM = FMBW3 
BL = BBW3 
XI = BX BS 
X2 = BXCS 
NH = 5 
GO TO 292 

C CHECK FOR INTERSECTION WITH BREAKWATER NO.C-D 
2 7 A FLM = FMBW2 
BL = BBW2 
XI = BXCS 
X2 = BXDS 
NH = 6 
GO TO 292 

C CHECK FOR DIFFRACTION ZONE OF END D. 

276 X 2 = BXDS + XDIST 
Y 2 = BYDS - YDIST 
XI = BXDS 
Y 1 = BYDS 
NH = 7 

USE SLOPE I NSER SECT I ON (M*X + B) TO SEE IF RAY AND 
BREAKWATER COMPONENT INTERSECT. 

290 IF ( ARS ( X 2-X 1 ) -0.00001) 2901,2902 *2902 
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2901 X 2 = X 2 + P.onoi 

2902 FLM = ( Y2-Y1 ) / ( X2-X1 ) 

PL = Y1 - FLM*X1 

292 IF ( ABS ( FLMR-F LM ) - 0.00001) 2°21 * 2922 *2922 

2921 FLMR = FL MR + 0.0001 

2922 XI = ( BL-BR > / ( FLMR-FLM) 

I F ( (X-XI ) * ( X I - X,N ) ) 294 *293*293 

293 I F ( ( X2-X I ) * ( X I-Xl ) ) 294*296*296 

2 q 4 GO TO ( 257*264,268 *274,274*276*92 ) *NH 
C 

C IF THEY DO INTERSECT CALCULATE POINT OF INTERSECTION VALUES* 
C INCLUDING RAY ANGLE* RAY SEPARATION FACTOR* 

C COEFFICIFNT OF REFRACTION AND SHOALING FACTOR. 

C 

296 Y I = ( X I-XN ) * ( Y-YN ) / ( X-XN ) + YN 

DFLS = ( ( X I -X N ) ** 2 + ( YI-YN ) **2 > **0. 5/S 
B1 = BNFWFI RN *B1 *B2 *DFLS ) 

A = BNEWF( AMI , AN, A, DFLS) 

AD = A * 180. /PI 
FKD = ( 1./B1 )**0. 5 
X = XI 
Y = Y I 
NN = 1 

CALL RSCAFK(NN) 

IF (NP-1) 2962*2962*2961 

2962 D = 1. 

GO TO 2963 
2961 NP = 3 

CALL RSBETA (NP,D*FKD) 

2963 WRITE (7*298 ) X I , Y I 
WRITE (7,297) 

WRITE (7,993) 

WRITE (7*998) N * X * Y , AD * C * B 1 * FKD » D * FK , S , DX Y 
298 FORMAT (41H INTERSECTION WITH BREAKWATER ZONE AT X= * 

1 F6.3 *4H Y= , F6 . 3 ) 

297 FORMAT ( 36H RAY VALUES AT INTERSECTION FOLLOW. ) 

GO TO 100 

92 GO TO (921*1°) »NDX Y 
CALCULATE BETA 

920 FKN = FK 

921 CALL RSBFTA ( NP , D » F KD ) 

GO TO 20 

END OF RAY CALCULATIONS. 

100 '/'RITE (7*994) 

994 FORMAT (13H RAY STOPPED ) 

RFTURN 
END 
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3. SUBROUTINE RSINTP 



Subroutine which calculates depth and partial derivatives of the bottom. 

SUBROUTINE NAME: RSINTP 

Variables used in Program 

COMMON variables are listed under MAIN 



DMMY ...... 

JX, JY . . . . 
DLX, DLY . . . 
IX , IY . c . . . 
NX, NY. . c . . 
IXNY , IYNY . . . 

D(NY) 

DD (NX) . . . . 
DD1 (NX) . . . . 
DD2 (NX) . . . . 



Dummy variables, no significance 
Integer X and X; not rounded. 

Decimal fraction of (X-JX) and (Y-JY) . 
(JX-2) and (JTY-2) . 

Counters for X and Y . 

X and Y grid points. 

Bottom depth at location DMAT (IXNY , IYNY). 
Bottom depth at (IXNX, Y). 

9d/3Y at (IXNX, Y) . 

9 2 d/8y 2 at (IXNX, Y) . 



Functions used in Program 

DEPTHF(A1, A2, A3, A4 , DELTA) 

Calculates data value at point of interpolation. Given, 
four equally spaced data values, Al, A2 , A3, and A4 , 
this function uses continuous parabolic interpolation, 
as described in the text, to calculate a data value at 
point DELTA. DELTA is that fraction of the distance from 
point A2 towards A3 at which we desire data values. 
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DEP1F (Al, A2 , A3, A4 , DELTA) 



Calculates slope of curve at point of interpolation 
using continuous parabolic interpolation* 

DEP2F (Al, A2, A3, A4 , DELTA) 

Calculates second derivative of curve at point of 
interpolation using continuous parabolic interpolation* 

Summary of Program 

Continuous parabolic interpolation is reduced to three functions; one 

to calculate a data value; one to calculate a slope; and one to calculate 

a second derivative. All values are at any abcissa value between the 

center two of four equally space data points* 

With (X,Y) located within the center grid of a 4 x 4 data array, 

these functions are used by 4 interpolating parabolas, parallel to the Y 

axis, to define 4 points on a plane parallel to the X axis and passing 

through (X, Y) . Using these 4 new points, and their derivatives when 

parabola 

necessary, the same functions are then used by an interpolating/parallel 
to the X axis, and passing through (X, Y) to completely define (X, Y) 
and all desired derivatives of a surface through (X, Y) . For a detailed 
explanation of this process see Section II-E (C) of text. 
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MACRO FLOW CHART OF SUBROUTINE RSINTP 



Enter from RSRFTN 
or RSCAFK 






MisCc Bookkeeping 
and setting constants 



Calculate values and derivatives 
for 4 tf Y intercepts' 1 



1 


l 


Calculate values 
for (X, Y) using 
values 


and derivatives 
4 "Y intercept” 



i 

^ Return 
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non nnnn nnnr> 



SUBROUTINE RSINTP ( DMMY ) 

C 

C SUBROUTINE CALLED BY RSRFTN OR RSCAFK WHICH WILL USE 
C 3-DIMENSIONAL PARABOLIC INTERPOLATION TO CALCULATE 

C DEPTH. BOTTOM SLOPES, AND OTHER PARTIAL DERIVATIVES 

C OF BOTTOM. 

C 

DIMENSION DM AT ( 50 , 40 ) » D ( 4 ) » DD ( A ) * DD 1 ( A ) >DD2 ( 4 ) 

COMMON DMAT»IXN»IYN»PI»G»GRID,X»Y,DXY,D1X,D1Y,D2X,D2Y,D2XY, 

1 A , AD, T »FK»BN»B1 »B2»C»PK»PPK»CPPK,CN»DDDC*DCDX»DCDY, 

2 S »NBW* PXAS *BYAS*BXBS*BYBS.BXCS *BYCS t BXDS *BYDS »WLD 

BOOKKEEPING AND SETTING OF CONSTANTS. 

FUNCTIONS FOLLOW 

DEPTHF ( A 1 » A2 » A3 ,A4»DELTA ) = A2 + ( DELT A / 2 . ) * < A3- A 1 ) - ( DEL T A** 2 /2 . ) 
1 * ( A4-4. *A3+5 .*A2— 2.*A1 )+(DELTA**3/2.>*(A4-3.*A3+3.*A2-Al) 

DEP1F(A1»A2,A3»A4,DELTA)= ( A3-A1 )/2. -( DELTA/2. ) *< A4-5. *A3+7. *A2 
1 -3. *41) + ( DE L T A**2 ) * ( A4-3 . *A 3 + 3 . * A2- A1 ) 

DEP2FI A.l ,A2 »A3 ,A4 , DELTA ) =A3-2 .*A2+A1 +DELT A* ( A4-3 . *A3 + 3 . *A2-A1 ) 
JY = Y 

DLY = Y- ( FLOAT ( JY ) ) 

JX = X 

DLX = X - ( FLOAT ( JX ) ) 

IX = JX -2 
IY = JY - 2 

INTERPOLATE 4 TIMES PARALLEL TO Y AXIS FIRST 

CALCULATE DEPTHS AND DERIVATIVES FOR Y-INTERCEPTS OF (X,Y) 

DO 200 NX = 1 >4 
I XNX = IX+NX 
DO 180 NYU, 4 
IYNY = IY + NY 
180 D(NY) = DMAT ( I XNX , I YNY ) 

DD (NX) = DFPT HF ( D ( 1 ) » D(2)» D(3>» 0(4), DLY) 

DD 1 ( NX ) = DEPIP ( D ( 1 ) » D ( 2 ) , D(3)» D(4)» DLY) 

200 DD2 ( NX ) = DEP2F (DU)* D(2), D(3>* D(4), DLY) 

USE RESULTS OF LAST 4 INTERPS TO CALC VALUES AT (X,Y) 

DXY = DEPTHF( DD( 1 ) ,DD ( 2 ) ,DD ( 3 ) ,DD( 4 ) ,DLX ) 

D1Y = DEPTHF ( DD1 ( 1 ) * DD1 ( 2 ) * DD1 ( 3 ) , DD 1 ( 4 ) , DLX ) 

D2Y = DPPTHF ( DD2 ( 1 ) ,DD2 ( 2 ) »DD2 ( 3 ) ,DD2 ( 4 ) ,DLX ) 

D 1 X = DEP1F (DD( 1 ) * DD ( 2 ) * DD ( 3 ) » DD ( 4 ) , DLX ) 

D2X = DEP2F ( D D ( 1 ) »DD(2) * DD ( 3 ) » DD ( 4 ) ,DLX) 

D2XY = DFP1 F ( DD1 ( 1 ) ,DD1 ( 2) ,DD1 ( 3 ) ,DD1 < 4 ) ,DLX ) 

400 RETURN 
END 
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4. SUBROUTINE RSBETA 



Subroutine which calculates the ray separation factor (BETA) , the 
coefficient of refraction, , and the shoaling coefficient, D, for a 
ray. 



SUBROUTINE NAME: RSBETA 

Variables used in Program : 

COMMON Variables are Used under MAIN. 



CG . . . . 
CGN . . . . 

D 

D2CDX2 . . 
D2CDY2 . . 
D2CDXY . . 
D2DDC2 . . 
FKD . . . . 
FPD . . . . 
PP, QQ • . 



Wave group velocity. 

Deep water wave group velocity. 

Shoaling coefficient; D in text, 

9^C/3X^ in text, 

3^C/3Y^ in text, 

3^C/3X3Y in text c 
2 2 

d d/dC in text. 

Coefficient of Refraction; in text. 

Convenient summation. 

p and q of equation for B2 shown under summary of program. 



Summary of Program 

After the ray trace subroutine, RSRFTN, has determined an (X, Y) 
coordinate, RSBETA is called to use a "Finite Difference" method, (See 
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text, Section H-E (d)) , for projecting from the two most recent Beta 
values to the next succeeding value. 

The program requires 2 initial Beta values to start, one at the 
present (X, Y) point, Bl, and one at the preceeding (XN, YN) point, BN, 
All points must be equally spaced. 

Each time called, RSBETA increments itself from the prior BN, Bl 
values and calculates the succeeding B2 using 



B2 



4-2 q s 2 

1 2 + pS 



•) Bl + (■ 



p S - 2 
2 + p S 



■) BN 



It must be kept in mind that this method is projecting an 
estimated Beta value one step ahead of ray tracing. One return from 
RSBETA, Bl is the present (X, Y) Beta value, 

In addition, the program calculates the Shoaling Coefficient at a 
point by using 



^ , CGNs 

D ’'cT 1 



1/2 



where 



CG = 0 c 5 * C * (1 + 



4 7T d/L v 
sinh (4 TT d/L) 



and CGN is the deep water value of CG. 

It also uses the Beta value to calculate the coefficient of 
refraction at a point. 



- 1 , 



K, = (Beta) 
d 



In program notation this is: 

FKD * t,Bl)” 1/2 
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MACRO-FLOW CHART OF SUBROUTINE RSBETA 
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non nnnn non n n n non non 



SUBROUTINE RS BET A ( NP » D » FKD ) 

C SUBROUTINE TO CALCULATE RAY SEPARATION FACTOR (BETA) FOR 
DIMENSION DM A T ( 50 , 40 ) 

COMMON DMAT ♦ I XN» I YN »PI »G»GRI D »X *Y»DXY»D1X »D1Y »D2X »D2Y» D2XY * 

1 A , AD * T ,FK,BN,B1 , B2 , C » PK » PPK , CPPK * CN , DDDC » DC DX , DCD Y , 

2 S»NB'.V*BXAS»BYAS»BXB^*BYRS»BXCS *BYCiJ »BXoS»BYDS»WLD 

INCREMENT BETA TO PRESENT POINT 

GO TO ( 700,700 ,705 ) »NP 
700 BN = B1 
B 1 = B2 

STILL IN RELATIVELY DEEP WATER. ..BETA = 1.000000 

IF (ABS(CPPK-l.O)-O. 0000001) 704,702,702 

704 B2 = 1. 

GO TO 710 

CALCULATE NECESSARY DERIVATIVES 

702 D2DDC2 = P K* ( 2 • * ( PPK/ ( l.+CPPK) +PPK/ ( l.-CPPK) ) + C*(PPK**2/ 

1 (l.-CPPK)**2 - PPK**2/ ( 1 .+CPPK ) **2 ) ) 

D2CDX2 = D2X/DDDC - ( D2DDC2*D 1X**2 ) /DDDC**3 
D2CDY2 = D2Y/DDDC - ( D2DDC2*D1Y**2 ) /DDDC**3 
D2CDXY = D2XY/DDDC - ( D2DDC2*D] X*D1 Y ) /DDDC**3 

SOLVE BETA EQUATION 

PP = - ( DCDX *COS ( A ) + DCDY*SIN( A) )/C 

00 = (D2CDX2*<SIN( A)**2)- ( D2CDX Y*2 . *S I N ( A ) *COS ( A ) ) 

1 + ( D2CDY2* ( COS ( A ) ) **2 ) ) /C 

B2 = ( (PP*S-2. >/(PP*S+2. ) )*BN + ( ( 4 . -2 . *QQ*S**2 ) / ( PP*S+2 . ) ) *B 1 
ON RETURN PRINT OUT B1 AS BETA FOR POINTS X,Y 

CALULATE COEFFICIENT OF REFRACTION. (FKD) 

FKD = ( 1 . /PI ) **0.5 

CALCULATE SHOALING COEFFICIENT. (D) 

705 FPD = 4.*PI*DXY/( C*T ) 

CG = 0. 5*C*( 1 + ( FPD/ ( 0.5*( EXP ( FPD ) -EXP ( -FPD ) ) ) ) ) 

GO TO ( 706.708,708 ) »NP 

706 CGN = CG 
D = 1. 

NP = 2 
GO TO 710 

708 D = ( CGN/CG ) **0. F 

710 RETURN 
END 
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5- SUBROUTINE RSCAFK 



Subroutine which calculates wave velocity and ray curvature at any 
given (X, Y) point. 

SUBROUTINE NAME: RSCAFK 

Variables Used in Program 

COMMON variables are listed under MAIN. 



CL Wave speed used to check convergence of iterations. 

LL ...... Index; no significance. 

NN o Index to indicate program status. If NN=1 , call 

RSINTP to determine DXY, then only calculate C. If 
NN=2, use present DXY and calculate both C and FK. 



Summary of Program 

When called with NN=1, RSCAFK calls RSINTP to determine water depth, 
then determines wave speed using an iteration process to solve: 



C 



= 

2tt 



tanh 



27Td 

L 



The program then returns. 

When called with NN=2, depth data have already been determined. 
Using these data, it determines wave speed then goes on to determine ray 
curvature using 



1 , . A , 9C n A ,3Cv N 

FK = — (sin A (g^) - cos A C ) 
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The program then returns. 

As a safeguard, a maximum of 80 iterations are allowed for C to 
reach successive values within 0^0005 of each other; otherwise a "Wave 
Speed Approximation" message is printed and the program continues using 
an average C value. 
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MACRO FLOW CHAT OF SUBROUTINE: RSCAFK 
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SUBROUTINE RSCAFK ( NN > 



C 

C SUBROUTINE WHICH CALCULATES WAVE VELOCITY AND RAY CURVATURE 
C AT ANY GIVEN (X»Y) POINT 

C USUALLY CALLED BY RSRFTN DURING RAY TRACE COMPUTATIONS* BUT 
C MAY BE CALLED BY MAIN JUST TO CALCULATE WAVE VELOCITY 

C AT ANY SPECIFIED LOCATION. 

C 

DIMENSION DMAT (50,40) 

COMMON DMAT , IXN,IYN»PI*G»GRID»X»Y,DXY»D1X*D1Y»D2X»D2Y»D2XY» 

1 a*ad*t,fk,bn*bi »82*c»pk»ppk*cppk»cn»dddc*dcdx»dcdy» 

2 S»NBW»BXAS»BYAS»BXBS »BYrS»BXcS»BYCS»BXDS»BYDS»WLD 

IF CALLED BY MAIN, NEED TO HAVE DEPTH CALCULATED FIRST. 

GO TO ( 5 7 , 55 > »NN 
53 CALL RSINTP ( DMMY ) 

CN = G * T / ( 2 .*P I ) 

C = CN 

CALCULATE WAVE SPEED, C. 

55 CL = C 
DO 56 LL = 1,80 

C = CN * TANH ( ( 2. *P I*DXY ) / ( T*CL ) ) 

IF (ABS(CL-C) - 0.0005) 58,58,56 

56 CL = ( C + CL ) /2. 

WRITE (7,999) 

999 FORMAT (26H WAVF SPFFD APPROXIMATED. ) 

IF CALLED BY MAIN. ..GO BACK TO MAIN 

58 GO TO ( 59,57 ) »NN 
CALCULATE CURVATURE OF RAY (= FK ) 

57 CP PK = C * PPK 

DDDC = (PK)*( (2.*CPPK)/( l.-CPPK**2 )+ALOG (l.+CPPK) 

1 - ALOG (l.-CPPK) ) 

DCDX = D1X/DDDC 
DCDY = DIY/DDDC 

FK = ( 1 . /C )*( ( SIN ( A) ) *DCDX - ( COS ( A ) ) *DCD Y ) 

59 RETURN 
END 
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6. SUBROUTINE RSDIFF 



Subroutine to calculate Coefficient of Diffraction and Direction of 
Wave Front Travel at end of semi-infinite breakwaters or at a breakwater 
gap. 

SUBROUTINE NAME: RSDIFF 



Variables used by Program : 

COMMON variables are listed under MAIN* 



AA , BB i « • • 
AA1, BB1 . . * 
ABW 

ABW1 

AF ..«.«« 

Al, A2 

BIX, B2X . . . 

BXA, thru 
BYD 



A, B in text® 

Used to save ’’half" solution across breakwater gap. 

Bearing of breakwater or gap, in radians, relative to 
X axis . 

ABW value used for direction checks. 

Bearing from breakwater tip or center of gap to (X, Y) 
and relative to X axis. 

Interim summations. 



For diffraction at one end of breakwater, the desired 
end is assigned coordinates, (BXB, BYB) ; the other end 
is assigned coordinates (BXA, BYA) . For diffraction at 
a gap, imaginary coordinates BXA through BYD are assigned 
based on an imaginary gap normal to the wave ray (see 
text) , 
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cc ..... 




C in text, part of Fresnel lntegral c 


CDIF . . . . 




Coefficient of diffraction,. 


CTM, CTP . . 


• 


cos (TH-THN) , cos (TH^THN) respectively. 


DADTH, DADR . 


( 


3A/30 and 3A/3r in text. 


DADX, DADZ . 


o 


3A/3X, 3A/3Z in text. 


DADZ1, DAD XL, 






DBDZ1, DBDX1 


• 


Used to save "half" solution across breakwater gap. 


DAI, DA2 , 






DB1 , DB2 . . 


c 


Interim summations 0 


DBDTH , DBDD . 


e 


3B/30 and 3B/3r in text. 


DBDX, DBDZ . 


e 


3B/3X, 3B/3Z in text. 


DFM, DFP , 






DU1, etc. . . 


a 


Appropriate derivatives used to make passes through 






similar equations without rewriting equations. 


DFMDR, DFP DR 


• 


3FH/3r, 3FP/3r in text. 


DFMDTH, 






DEPDTH . . . 


• 


3FM/39, 3FP/30 in text. 


DG ..... 




Gap length. 


DIR 




Angle of wave travel in radians. 


DIRD .... 




Angle of wave travel in degrees. 


DSIG2R . . . 


• 


3a f /3r; also 3a/3r; see Note. 



Note: Same calculations are needed with both a and o ! . 

A first pass is made through calculations (M=l) using 
a values although the variable name indicates b T values. 
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DSIG2T . . 

DT . . . . 

DU2DR, 
DW2DR . . . 

DU2DTH, 
DW2DTH . . 

DX, DY . . 

DXT , DYT . 

FACFP , 

FBCFB . . . 

FACFP X, 
FBCFPX . . 

FACFPZ , 
FBCFPZ . . 

FFl , FF2 . 

FLTK . . . 

FM, FP . . 



Upon completion of calculations o, results are saved, o' 
values are inserted, and a second pass is made (M=2) . 
First pass, with o, calculates U1 , W1 and their deriva- 
tives; second pass, with o', calculates U2, W2 and their 
derivatives . 

. 9a’/96; also 9a/90; see Note following DSIG2R. 

. Half width of imaginary gap normal to wave ray. 

. 3U2/9r , 9W2/9r or 3Ul/9r, 9Wl/3r. See note under 

DSIG2R. 

. 9U2/90, 3W2/30 or 9U1/90, Wl/30. 

. Various X and Y differences used locally in various parts 
of program; no significance. 

. Local subtotal; no significance. 

. cos (FP) , sin (FP) in text. 

. 3 (cos (FP) ) /9X, 9 (sin(FP) ) /3X in text. 

. 3 (cos (FP) ) /3Z , 9 (sin(FP) ) /3Z in text. 

. Interim summations. 

. K in text; 2 *PI/(C*T). 

. FM, FP in text; kr cos (TH-THN) , kr cos (TN+THM) 
respectively. 
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FN 



n in text, index for Fresnel summations. 



2 

FT ...... t in text; Tfo 

Fl, F2 .... Variables from MAIN; coordinate locations (X, Y) c 

IWD2 . . . . . Alphameric variable from MAIN; the end of breakwater 
which is being diffracted, A, B, C, or D. 

LL . . o c . o Index passed from MAIN which indicates desired action by 
subroutine. If LL~1, 2, 3, or 4, diffraction is at A, B, 
C, or D end of breakwater respectively. If LL=5 , diffrac- 
tion is at Gap 0 If LL=6 or 7, perform diffraction for 
given coordinates. 



M . Index; see Note following DSIG2R. 

N Index (1, 2, or 3) used to return to correct location in 



program after calculating a bearing. 

NGAP ..... NGAP=1 for diffraction at the end of a breakwater. 

NGAP=2 for gap diffraction with higher values being 
assigned during the process of gap solution in order to 
keep track of the status of calculations. 



NGAPS Original value of NGAP to be restored at end of program. 

NN NN=1 if imaginary gap is being established at start of 

program; NN=2 at end of program when gap is being restored 
for future use. 



NQ . Index to indicate whether angle AP is greater than 

angle A (NQ=2) , or less (NQ=1) . 

R r in text; distance from (X, Y) to either a breakwater 

tip or center of gap c 
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SI 


o or o 1 . 


SI4 


4 

O o 


SIGI 


o in text. 


SIG2 


o* in text. 


SS 


S in textc Fart of Fresnel Integral. 


TERM 


Terms of summation series 0 


TH 


Angle from breakwater to (X, Y) as measured at breakwater 
tip. 


THN 


Angle from breakwater to ray of incident wave as 
measured at tip 0 


THD, THND 


TH and THN in degrees. 


THX 


0.5 (0+0q) or 0.5 (0-0q) ; no significance. 


TH11 


Used to save "half" solution across breakwater gap. 



Ul, U2, 



Wl, W2 


Same as texto 


XC, YC 


Coordinates of center of gap. 


ZZ 


0.5 (2*R/(C*T)); no significance. 


Z2 


Convenient summation. Used in sin and cos approximations 
of Fresnel Integrals. 
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Summary of Program 



Prior to starting diffraction calculations, incident wave data must 
be input to the program as previously described in the section on MAIN 
commands . 

In order to solve a diffraction problem, RSDIFF is first called by 
MAIN after a DIFFRACT Command has been read. Depending on whether told to 
diffract breakwater end A, B, C, D, or Gap, index LL is set to 1-5 
respectively by MAIN. Then RSDIFF is called to assign breakwater co- 
ordinates (BXB , BYB) to the end of a breakwater at which diffraction is to 
take place and assigns BXA, BYA to the other end. BXAS through BYDS are 
never disturbed and are kept for later use if desired. The main portion 
of RSDIFF is such that it always computes diffraction about the B end of 
an A-B breakwater. When called to diffract a gap, the program first 
determines an imaginary gap, (see Text Section III, B) which is normal 
to the incident wave ray and sets index NGAP=2 to indicate a gap problem. 
Control then returns to the main program. 

When later called by MAIN after a DIFFR COORD command has been read, 
RSDIFF calculates the coefficient of diffraction and direction of wave 
travel at (X, Y) and about the previously designated tip or gap. 

The main portion of the RSDIFF program is arranged to calculate 
values needed to solve a standard semi- inf inite breakwater problem. 
Breakwater gap solutions can be built out of semi- inf inite solutions. To 
solve for a gap, it is therefore necessary to make multiple passes through 
the body of the program with "semi-infinite breakwater" portions of a gap 
then add the results for the final desired values. 

Equations for this program are explained in detail in the text of 
this report. Due to their length and complexity, they are not repeated 
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in the summary which follows; but to the maximum extent possible, similar 
notation is used in the programs and in the text to facilitate cross 
reference. 

Using text notation TH, THN, the program first determines necessary 
angles, and R. These are then used to find FM, FP and their theta and R 
derivatives . 

Next SIGMA 1 and SIGMA 2 are derived and used in Fresnel integrals to 
determine C and S* These are then used to give Ul, W1 , U2, W2 and the 
derivatives of each with respect to both Theta and R. 

At this point, one has all the input needed to find A, B, and their 
Theta and R derivatives. 

Continuing the building process, the derivatives of A and B are used 
to find 3A/8X, 8A/3Z, 8B/8X, and 8B/8Z C 

If the solution is only desired for a semi-infinite breakwater, 
these values are used directly to find the coefficient of diffraction and 
the wave direction of travel. 

If a gap solution is required, A, B, and their derivatives are saved 
as a partial solution and coordinates (BXB, BYfe) are established to coin- 
cide with the "center” of the breakwater gap. It is then necessary to go 
back and repeat those steps necessary to solve for FP and its derivatives. 
This information is used to determine exponential derivatives needed for 
the solution. 

Coordinates (BXB, BYB) , (BXA, BYA) are then changed to coincide 
with the C-D breakwater. All steps are repeated as required to find A, 

B, and all their X, Z derivatives again. 

It is then possible to add all partial solutions for the gap to give 
totals for A, B, and their derivatives. These are used to find the 
coefficient of diffraction and the wave direction of travel. Upon 
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i 














completion of a gap solution, coordinates of the imaginary breakwater are 
restored so that further calculations can be made if desired. 
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MACRO-FLOW CHART OF RSDIFF 



c 



Enter from MAIN 






A,B,C,D 



Go to location specified 
by second character of command 



f 



GAP 



Set up (BXB ,BYB) 
and (BXA,BYA) to 
perforin diffrac- 
tion off "B" end 
of breakwater. 



R 



I 



i 



Establish imaginary gap 
normal to wave ray 
Set up diffraction off 
"B" end. 



pZ 

JL 



Set up coords toj 
do calc, off 
center of gap 



Set X,Y | 



Determine 
AP , ABW , R 



Determine 

TH,THN 



Calc. FM, FP 
and their 
derivatives 




Save values 
Set up to 
do calc. off 
"C" endo 



Restore Gap 
if used 



\ Output 7 
\ Results j 



Calc . SIGMA2 
and its 
derivatives 



Calc. SIGMA1 and 
its derivatives 

---L 



i 



Evaluate Fresnel Integrals 
to get C and S 



Save Ul, 
Wl, etc. 



i 



Calc. U2, W2 and 
their derivatives 



Calc. K T and Wave 
Dir. of Travel 








Save old 
values 









Calc. A, B, and their 
derivatives (Two passes 
through same equations) 



i 



Calc. DADX, DBDX, 
DADZ, DBDZ 
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SUBROUTINE RSD I FF ( LL , I WD2 * F 1 * F2 ) 

C PROGRAM TO CALCULATE COEFFICIENT OF DIFFRACTION AND DIRECTION OF 
C WAVE (FRONT) TRAVEL AT END OF BREAKWATERS OR AT A 

C BREAKWATER GAP. 

DIMENSION DMAT ( 50 ♦ 40 ) 

COMMON DMAT* IXN, IYN»PI , G , GR I D , X , Y * DXY * D IX , D1 Y » D2X , D2 Y * D2 XY » 

1 A»AD»T»FK»BN,B1»B2»C»PK»PPK»CPPK»CN»DDDC»DCDX»DCDY» 

2 s,nbw*bxas,byas*bxbs»bybs*bxcs*bycs»bxds»byds*wld 

C 

C GO TO LOCATION SPECIFIED BY COMMAND . 

C 

GO TO (5001. 5002*8003*5004, 5005, 5006*5006) ,LL 
C 

C BOOKKEEPING TO SET UP DATA FOR CALCULATIONS 
C PROGRAM CALCULATES VALUES OFF B END OF BREAKWATER 
C FOR SOLVING GAP, THREE PASSES WILL BE NEEDED TO SOLVE 
C FOR THREE IMAGINARY ENDS . 

C 

C DIFFRACTION AT A END OF BREAKWATER 

5001 BXB = BX AS 
BX A = BXBS 
BYB = BYAS 
BYA = BYBS 
NGAP = 1 
GO TO 5008 

C DIFFRACTION AT B END OF BREAKWATER 

5002 BXB = BXBS 
BX A = BXAS 
BYB = BYBS 
BYA = BYAS 
NGAP = 1 
GO TO 5008 

C DIFFRACTION AT C END OF BREAKWATER 

5003 BXB = BXCS 
BYB = BYCS 
BX A = BXDS 
BYA = BYDS 
NOAP = 1 
GO TO 5008 

C DIFFRACTION AT D END OF BREAKWATER 

5004 BXB = BXDS 
BYB = BYDS 
BXA = BXCS 
BYA = BYCS 
NGAP = 1 

5008 WRITE (7,5099) IWD2 

5099 FORMAT (30H DIFFRACTION CALCULATIONS FOR , Al, 

1 19H END OF BREAKWATER. ) 

RETURN 

C 

C DIFFRACTION IN BREAKWATER GAP 



c establish an imaginary gap. 

c 

c 



5005 

5009 



5007 



50 L 5 
5098 

5097 

5999 

5016 



NGAP 
NN = 

DX = 

DY = 

DG = 

XC = 

YC = 

N = 3 
GO TO 503 
DT = DG/2. 



= 2 
1 

BXCS - 
BYCS - 
( DX**2 
( BXBS h 
( BYBS H 



BXBS 

BYBS 

+ DY**2)**0.5 
BXCS)/2. 

BYCS ) /2 ♦ 



DXT = 
DYT = 
BXB = 
BX A = 
BXC = 
BXD = 
BYB = 
BYA = 
BYC = 
BYD = 
GO TO 
WRITE 
FORMAT 



* SIN(A-ABW) 
SIN(A) 

COS ( A ) 

DXT 

3.*DXT 
DXT 

3 • *DXT 
DYT 

3.*DYT 
DYT 

3.*DYT 
< 5015*5016) *NN 
(7,5098) 

(46 H DIFFRACTION 



DT 

DT 

XC 

XC 

XC 

XC 

YC 

YC 

YC 

YC 



# 

* 



+ 

+ 

+ 

+ 



CALCULATIONS FOR BREAKWATER GAP. 



WRITE (7,5097) 

FORMAT (52H BASED ON IMAGINARY GAP ASSUMED MORMAL TO WAVE 
WRITE (7,5999) BXB , BYB , BXC , BYC 

FORMAT ( 2 1H IMAGINARY GAP FROM , F6 . 3 , 1H , , F6 . 3 ♦ 

4H TO ♦ F6 • 3 * 1H » , F6 . 3 ) 

RETURN 



C 

C SET COORD POINTS 
C 



5006 X = FI 
Y = F2 

NGAPS = NGAP 
C 

C DETERMINE BEARING AND DISTANCE FROM BWB TO X,Y (AP,R) 
C 

500 N = 1 

DY = Y - BYB 
DX = X- BXB 

R = ( DX**2 + DY**2 ) **0 • 5 * GRID + 0.0000001 

503 IF ( ABS(DX)-0. 00001 ) 5011,5012,5012 

5011 ABW = PI/2. 

GO TO 5014 

5012 ABW = AT AN ( DY/DX ) 

C QUADRANT CHECK 

5014 IF (DY*DX) 504,504,506 



) 

RAY. ) 
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504 IF(DX) 507,505,508 

505 ABW = PI/2. 

506 IF(DV) 507,509,500 

507 A8W = ABW + PI 
GO TO 509 

508 ABW = ABW + 2.*PI 

IF ( ARW-2 • *P I ) 509,502.502 
502 ABW = ABW - 2.*PI 

509 GO TO (511,513,5007) ,N 

DETERMINE BREAKWATER BEARING (ABW) 

511 AP = ABW 

DY = BYA - BYB 
DX = BXA - BXB 
N = 2 
GO TO 503 

DETERMINE ANGLE FROM BW TO RAY(THN) 

513 THN = ABS ( A-ABW ) 

DETERMINE ANGLE FROM BW TO X»Y (TH) 

TH = ABS(AP-ABW) 

IF (TH-PI) 516,516,515 

515 TH = 2.*PI - TH 

516 IF (THN-PI) 518,518,517 

517 THN = 2 • *P I - THN 

SET MISC SUMMATIONS TO BE USED LATER 

518 FLTK = 2.*PI/(C*T) 

ZZ = (FLTK*R/PI ) **0 • 5 
CTM = COS (TH-THN) 

C TP = COS (TH+THN) 

CALC FM, FP, AND THEIR DERIVATIVES. 

FM = FLTK*R*CTM 

FP = FLTK*R*CTP 

DFMDR = FLTK*CTM 

DFPDR = FLTK*CTP 

DFMDTH = -FLTK*R*SIN (TH-THN) 

DFPDTH = -FLTK*R*S IN (TH+THN) 

NGAP = 1....NO BW,,, NGAP = 2. ...DOING LEFT SIDE 
NGAP = 3 ....DOING CENTER, ,♦♦ NGAP = 4 ....DOING RIGHT SIDE 
GO TO ( 519,519,577,519 ), NGAP 

USE SIGMA VALUES TO CALC U, W, AND THEIR DERIVATIVES 
THROUGH THE USE OF FRESNEL INTEGRALS. 

CALC SIG1 VALUES F I RST . . . . THEN SIG2 VALUES (USE LOOP) 
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519 SIG2 = 2 • *ZZ*S I N ( 0 • 5* ( T H-THN ) ) 

THX = 0 . 5 * ( TH-THN ) 

IF (SIG2) 5191 *5191*5192 

5191 NO = 1 

GO TO 5193 

5192 NO = 2 

5193 M = 1 

510 FT = PI*SIG2**2/2. 

IF (ABS(THX) - 0.000001) 5 1 0 1 * 5 1 02 » 51 02 

5101 DSIG2T = ZZ 

GO TO (5104.5105 ) .M 
5105 DSIG2T = -ZZ 
GO TO 5104 

5102 DSIG2T = SIG2*COS(THX)/(2.*SIN(THX) ) 

5104 SI = ABSISIG2) 

EVALUATE FRESNEL INTEGRALS ( CC AND SS ) 

IF (SI-3.) 512.540 »540 
IF SI L.T. 3. USE SERIES SOLUTION 
512 SI = Sl*( (PI/2. )**0.5) 

S 1 4 = SI**4 
FN = 1. 

TERM = SI 
CC = SI 

514 FN = FN+1 • 

TERM = - TERM*Sl4*(4.*FN-7. ) / ( 16 .*FN**3-52 . *FN**2+54.*FN-18 . ) 
CC = CC+TERM 

IF (ABS (TERM) - 0.0001) 524.514.514 
524 CC = CC*( (2./PI )**0.5) 

FN = 1. 

SS = (SI**3)/3. 

TERM = SS 
530 FN = FN+1. 

TERM = -TERM*SI4*(4.*FN-5. )/( 16 . *FN** 3-28 . *FN**2+14 .*FN-2 . ) 

SS = SS+TERM 

IF (ABS (TERM)-O.OOOl) 538.530.530 
538 SS = SS*< (2. /PI )**0.5) 

GO TO 545 

IF SI G.T. 3. USE SIN/COS APPROXIMATIONS 
540 Z2 = (Pl/2.)*(SI**2-2./<PI*SI)**2) 

CC = 0.5 + SIN ( Z2 )/(PI*SI ) 

SS = 0.5 - COS (Z2)/(PI*SI ) 

CALC U2.W2.AND ALL THEIR DERIVATIVES. 

545 U2 = 0.5*( l.-SS-CC) 

W2 = 0 • 5* ( SS-CC ) 

IF ( ABS( SIG2 ) -0.000001 ) 5451.5452.5452 
5452 DSIG2R = ( S 162/ ( 2 . *R ) ) * ( S I G2 / ( ABS ( S I G2 ) ) ) 

DSIG2T = DSI62T * ( S I G2 / ( ABS ( S I G2 ) ) ) 

^\^ 2 > - 
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GO TO 5453 

5451 DSIG2R = SIG2/(2.*R) 

5453 DU2DR = -0 . 5* ( S I N ( FT ) *DS I G2R + COS ( F T ) *DS I G2 R ) 

DW2DR =0.5*(SIN ( FT ) *DS I G2R-C0S (FT)*DSIG2R) 

DU2DTH = -0 • 5* ( S I N (FT)*DSIG2T + COS <FT)*DSIG2T) 

DW2DTH = 0 • 5* ( S I N (FT)*DSIG2T - COS <FT)*DSIG2T) 

ON FIRST PASS RESET VALUES TO SHOW U1»W1,AND THEIR DERIVATIVES 
SET SIGMA2 AND GO BACK TO CALCULATE U2,W2,AND THEIR DERIVATIVES. 

GO TO ( 550,560) ,M 
550 U1 = U 2 
W1 = W 2 
DU1DR = DU2DR 
DW1DR = DW2DR 
DU1DTH = DU2DTH 
DW1DTH = DW2DTH 

SIG2 = -2 • *ZZ*S I N ( 0 » 5* ( TH+THN ) ) 

THX = 0 . 5* ( TH+THN ) 

M = 2 
GO TO 510 

USE SAME BASIC A AND B EQUATIONS TO CALCULATE BOTH R AND THETA 
DERIVATIVES OF A AND B. 

FIRST DO R DEVIVATIVES. 

560 M= 1 

DFM = DFMDR 
DFP = DFPDR 
DU1 = DU1DR 
DU2 = DU2DR 
DW1 = DW1DR 
DW2 = DW2DR 

A 1 = U 1 * COS(FM) + W1 * SIN(FM) 

A2 = U2 *COS ( FP ) + W2 *SIN(FP) 

BIX = W1 *COS(FM) - U1 * SIN(FM) 

B2X = W2 * COS(FP) - U2 * S IN C FP ) 

565 DAI = DU 1 *COS(FM) - U 1*S I N ( FM ) *DFM + DW1*SIN(FM) 

1 + Wl*COS(FM)*DFM 

DA2 = DU2 * COS(FP) - U2*S I N ( FP ) *DFP + DW2 * SIN(FP) 

1 + W2*COS(FP)*DFP 

DB1 = DW1 * COS(FM) - W1*S IN ( FM ) *DFM - DU1*SIN(FM) 

1 - Ul*COS(FM)*DFM 

DB2 = DW2*COS ( FP ) - W2*S IN ( FP ) *DFP - DU2*SIN(FP) 

1 - U2*COS(FP)*DFP 

C WATCH OUT FOR PROPER QUADRANT 
GO TO (570,580) , NQ 
570 DADTH = DAI + DA2 
DBDTH = DB1 + DB2 
GO TO (572,575) *M 
572 AA = A1 + A2 
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BB = BIX + B2X 
GO TO 585 

580 IF(NGAP-l) 5801 *5801 *5811 

5801 DADTH = -S I N ( FM ) *DFM - DAI + DA2 

DBDTH = -COS ( FM ) *DFM - DB1 + DB2 

GO TO 5812 

5811 DADTH = DA2 -DAI 

DBDTH = DB2 -DB1 

5812 GO TO ( 582*575 ) »M 

582 IF(NGAP-l) 5821*5821*5822 

5821 AA = COS(FM) - A1 + A2 
BB = -SIN(FM) - B1X+ B2X 
GO TO 585 

5822 AA = — A 1 + A2 
BB = -BIX + B2X 

SAVE R VALUES.. . .SET VARIABLES SO AS TO CALCULATE THETA 
DERIVATIVES AND DO IT OVER AGAIN. 

585 M = 2 

DADR = DADTH 
DBDR = DBDTH 
DFM = DFMDTH 
DFP = DFPDTH 
DU1 = DU1DTH 
DU2 = DU2DTH 
DW1 = DW1DTH 
DW2 = DW2DTH 
GO TO 565 

CALC DERIVATIVES OF A AND B WITH REGARDS TO X AND Z. 

575 DADX = DADR*SIN (t H ) + DADTH*COS (TH) / R 
DBDX = DBDR*S I N (TH) + DBDTH*COS (THJ/R 

WATCH OUT FOR QUADRANTS 

5753 GO TO ( 57 5 1 ♦ 5 75 2 ♦ 575 1 » 575 1 ) ♦ NGAP 

5752 DADZ = -DADR*COS ( TH ) + DADTH*S I N ( TH ) / R 
DBDZ = -DBDR*COS(TH) + DBDTH*S I N ( TH ) /R 
GO TO 575 A 

5751 DADZ = DADR * COS(TH) - DAD TH*S I N ( T H ) /R 
DBDZ = DBDR*COS ( T H ) -DBDTH*S I N (TH)/R 

IF THIS IS A GAP GP BACK AND DO AGAIN FOR CENTER AND RIGHT SIDE 

5754 GO TO ( 574 ♦ 576 * 577 » 578 ) *NGAP 

SAVE OLD VALUES FOR LEFT SIDE OF BW 

576 DADZ 1 = DADZ 
DBDZ 1 = DBDZ 
DADX 1 = DADX 

- \ 3 T => - 
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DBDXl = DBDX 
TH1 1 = TH 
AA1 = AA 
BB1 = BB 

C DO CENTER NEXT 
NGAP = 3 
BYA = BYB 
BX A = BXB 
BXB = (BXB+BXO/2. 

BYB = ( BYB+BYC ) / 2 • 

GO TO 500 

SAVE CENTER VALUES (ONLY A PARTIAL SET) 



577 FACFP = COS(FP) 

FBCFP = -SIN(FP) 

FF1 = DFPDR*SIN( TH ) + DFPDTH * COS(TH)/R 
FF 2 = DFPDR*COS(TH) - DFPDTH*S IN ( TH ) /R 
FACFPX = -SIN(FP)*FF1 
FBCFPX = COS ( FP ) *FF 1 
FACFPZ = SINCFP)*FF2 
FBCFPZ = -COS ( FP ) *FF2 
GO BACK AND DO RT SIDE OF BW 
BX A = BXD 
BYA = BYD 
BXB = BXC 
BYB = BYC 
NGAP = 4 
GO TO 500 



IF THIS WERE A GAP PROPLEM MAKE SUMMATIONS. 



DOES POINT LAY BETWEEN GAP 
MAKE BW SUMMATIONS 



TIPS OR BEHIND A BREAKWATER. 



578 

581 



IF ( ( TH1 1 + 
DADX = DADX 
= DBDX 
= DADZ 
= DBDZ 
AA + 

BB + 



579 



DBDX 
DADZ 
DBDZ 
AA = 
BB = 
DADX 
DADZ 
DBDX 
DBDZ 
AA = 
BB = 



A - PI )*<TH - 
+ FACFPX 
+ FBCFPX 
+ FACFPZ 
+ FBCFPZ 
FACFP 
FBCFP 



A)) 



= DADX1 
= DADZ 1 
= DBDXl 
= DBDZ1 



DADX 

DADZ 

DBDX 

DBDZ 



AA1 

BB1 



AA 

BB 



MADE IT WHEW 

NOW CALC DIR. OF TRAVEL AND COEFF. OF DIFFRACTION 
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574 IF ( ABS< AA*DBDZ-BB*DADZ )-0 .00001 ) 5741,5742 .5742 

5741 DIR = PI/2. 

GO TO 5744 

5742 DIR = ATAN ( ( AA*DBDX -BB*DADX ) / ( AA*DBDZ~BB*C ADZ ) ) 

C DIR IS WRT BW RETURN TO X,Y COORD SYSTEM 

C CORRECT FOR NEGATIVE ATANS IN SECOND QUADRANT. 

5734 IF (DIR) 5743,5744,5744 

5743 DIR = DIR + PI 

5744 ABW1 = ABW 

IF (ABS(A-ABWI)-PI ) 5748,5748,5749 
5749 ABW 1 = ABW1-2 • *P I 

5748 IF ( (2. *PI-(A-ABW1 ) )*(A-ABW1 ) >5746,5746,5745 

5745 DIR = DIR + ABW 
GO TO 5747 

5746 DIR = ABW - DIR 

5747 IF (DIR-2.*PI) 5731,5732,5732 
5732 DIR = DIR-2.*PI 

5731 DIRD = D I R* 180 . /P T 

CD I F = (AA**2 + BB**2)**0.5 

PRINT OUT RESULTS 

WRITE (7,5096) 

WRITE (7,5095) AD,T»C 
WRITE (7,5094) X,Y 
WRITE (7,5093) D I RD» CD I F 
NGAP = NGAPS 

RESTORE ARTIFICIAL GAP IF USED 
NN = 2 

IF (NGAP-2) 5016,5009,5009 
5096 FORMAT ( 30H WAVE CHARACTERISTICS ) 

5095 FORMAT (9H ANGLE ,F7.2, 18H DEG., PERIOD »F6.2» 

1 16H SEC., SPEED, F5.2, 9H FT/SEC. ) 

5094 FORMAT (14H AT POINT X = , F5.2, 5H Y= ♦ F5.2) 

5093 FORMAT (27H DIR OF TRAVEL IS »F7.2, 

1 29H DEG.,COEF OF DIFFRACTION IS , F7.4) 

RETURN 
END 
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7 . SUBROUTINE RSSHOR 



Subroutine which calculates shoreline location from a depth array. 

SUBROUTINE NAME: RSSHOR 

Variables used in Program : 

COMMON variables are listed under MAIN. 

D • « Fraction of distance between two grid points at which 

shoreline crossing has been interpolated. 

DX, DY « « . . D distance along X or Y grid line respectively. 

INC • • . * ® * Index used to keep track of which of 6 possible grid 

crossings is being checked; INC=1 during original search 
for a shoreline segment. 

IX, IY .... Coordinates at primary point of shoreline search; 

program looks for shoreline between (IX, IY) and 
(IX + 1, IY) or between (IX, IY) and (IX, IY + 1). 

1X2, IY2 . . . Coordinates at secondary grid point opposite (IX, IY) . 

Shoreline check occurs between two points. 

NCK(LL) .... A check list of coordinate locations at which grid 
crossings have already been found. Coordinates are 
stored as an integer (100 * IX + IY) ; with a plus value 
indicating grid crossing between (IX, IY) and (IX + 1, IY) , 
and a negative value indicating grid crossing was between 
(IX, IY) and (IX, IY 4* 1). 
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NN 

NNS 



NS 



NSL 



XD, YD 
XS, YS 



Do loop counter, no significance. 

Index; NNS^l if last grid crossing was along Y grid line; 
NNS=2 if last grid crossing was along X grid line. 

Index; NS*=+ 1 if searching between IX and IX +1; NS= -1 
if searching between IY and IY +1. 

Index of total number of grid crossing located and 
stored in NCK(--) * 

Coordinates of shoreline at a grid crossing point. 

First coordinate points along a shoreline segment; used 
to close loop if shoreline traces back to origin as on 
an island or enclosed bay c 



Summary of Program 

When called by MAIN, RSSHOR finds and traces shoreline segments 
across a DMAT array by using the fact that, only at a shoreline, the 
product of two adjacent depths (DMAT (IX, IY) ) * (DMAT (1X2, IY2)) is 
negative. 

Starting at (2,2), the product of adjacent grid depths is determined 
along grid lines parallel to the Y axis until an initial starting point 
of a shoreline segment is found. Once a shoreline is found, it is then 
traced to its ending, either at the matrix edge or back at its start such 
as at an island. 

Each time a shoreline grid crossing is found, the coordinates of 
(IX, IY) are stored as integer (100 * IX + IY) in NCK (NSL). It is given 
a negative value if the grid crossing were between (IX, IY) and 
(IX + 1, IY) ; it is given a positive value if the crossing were between 
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(IX, IY) and (IX, IY + 1) . Each possible crossing point is first checked 
against NCK to be sure that it is in fact a new crossing point, not a 
duplicate. 

Shoreline traces are made by checking the six adjacent grid bounds 
for a new crossing; when found, it is stored in NCK. A check is then made 
of its 6 adjacent grid bounds. The NCK listing eliminates any line backing 
up upon itself. 

Each time a shoreline is ended, a new search starts at (1,2) until 
all segments have been found. 

Output is a series listing of consecutive shoreline coordinates. 
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MACRO-FLOW CHART OF RSSHOR 
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SUBROUTINE RSSHOR ( DMMY ) 

C SUBROUTINE WHICH WILL LOCATE THE SHORELINE 
C BASIC APPROACH. . .PRODUCT OF ADJACENT DEPTHS IS NEGATIVE 
C WHEN CROSSING SHORELINE) 

C USE STRAIGHT LINE INTERPOLATION BETWEEN GRID POINTS OF APPOSITE 
C SIGN TO APPROXIMATE SHORELINE 

C USE NCK (XXX) AS CHECK LIST FOR COORDINATE LOCATIONS 
C 

C BOOKKEEPING 
C 

DIMENSION DMAT( 50,40 > ,NCK( 200 ) 

COMMON DMAT,IXN»IYN»PI » G » GR I D , X , Y , DXY , D IX , D 1 Y , D2 X ♦ D2 Y , D2 XY » 

1 A* AD»T»FK*BN*B1*B2 , C » PK , PPK »C PPK , CN , DDDC , DCDX , DCDY » 

2 S,NBW,BXAS»BYAS»BX8S,BYBS»BXCS,BYCS,BXDS,SYDS»WLD 
NSL = 1 

C ZERO OUT CHECK LIST 
DO 402 LL = 1,200 
402 NCK(LL) = 0 

SEARCH FOR INCREMENT OF SHORELINE. 

INC = 1 .. NOT IN PROCESS OF TRACING SHORELINE SEGMENT 



404 


INC 


= 1 




IX 


= 1 


406 


NS 


= 1 




1X2 


= IX 




IY 


= 2 


410 


I Y2 


= IY + 1 


EDGE 


CHECK 


420 


IF 


( ( IX-1 )*< 


422 


IF 


( ( IY-1 )*{ 


CHECK 


PRODUCT 


425 


IF 


( DMAT ( IX, 



I XN- IX)) 490,422,422 
IYN-IY) ) 490,425,425 



IY)*DMAT< 1X2, IY2) ) 430,426,411 



INCREMENT TO NEXT POINT 

5EEP INCREMENTING UNTIL ENTIRE GRID IS CHECKED... 

411 IF (INC-1) 460,412,460 

412 IY = IY+1 

IF (IYN-IY) 414,414,410 
414 IX = IX + 2 

IF (IXN-IX) 495,406,406 

ELIMINATE ANY ZERO DFPTHS 

426 IF ( DMAT ( I X , I Y ) ) 428,427,428 

427 DMA T ( I X , I Y ) = DMAT(IX.IY) + 0.0001 

- \42 - 



nnnnnn nnn nnnnn nnnn nno 



428 IF ( DM AT ( I X2 » I Y2 ) ) 425,429,425 

429 DMAT < 1X2 , I Y2 > = DMAT ( I X2 » I Y2 ) + 0.0001 
GO TO 425 

CHECK FILE FOR PREVIOUS LISTING OF POINT 



430 DO 445 NN = 1,NSL 

IF (NCK(NN) - NS* ( 100*1 X+IY) ) 
445 CONTINUE 



445.411 .445 



INTERPOLATE BETWEEN GRID POINTS 
WAS IT AN X OR Y CROSSING... 





D = 


DMAT ( IX 


♦ IY 


) / ( DMAT ( : 


IX , I Y 




IF 


(NS-1) 450, 


452,452 




450 


DX 


= D 










DY 


= 0. 










NNS 


= 2 










GO 


TO 455 








452 


DY 


= D 










DX 


= 0. 










NNS 


= 1 








STORE 


GRID PO 


INT 


IN FILE 




STORE 


+ VALUE 


IF 


ALONG Y 


GRID 


STORE 


- VALUE 


IF 


ALONG X 


GRID 


455 


NCK 


(NSL) = 


NS* 


( 100*IX + IY ) 




XD 


= FLOAT ( 


IX ) 


+ DX 






YD 


= FLOAT ( 


IY) 


+ DY 






IF 


( INC-1 ) 


457 


,456,457 




456 


XS 


= XD 










YS 


= YD 









- DMAT ( 1X2 ♦ IY2 ) ) 



496 



WRITE (7.496) 

FORMAT (5X.22H SHORELINE 



COORDINATES) 



PRINT CROSSING POINT 



457 WRITE (7.497) NSL.XD.YD 
497 FORMAT < I5.17X.2F8.3) 

NSL = NSL+1 

INC = 2 WHEN FOLLOWING SHORELINE 
INC = 2 

SEARCH FOR ADJOINING SHORELINE GRID CROSSING 



LOOK 


IN 


6 POSSIBLE DIRECTIONS 


460 


GO 


TO 


( 459,470 ) , NNS 


459 


GO 


TO 


(462,462,463,464,465*466,467,468) , INC 


462 


IX 


= 


IX+1 
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n n n 



INC = 3 
461 1X2 = IX 
NS = 1 
GO TO 420 

463 IX = IX-2 

IF (IX) 4652*4652 *4631 
4631 INC = 4 

GO TO 461 

464 INC = 5 
IY2 = IY 

469 1X2 = IX + 1 
458 NS = -1 

GO TO 420 
4652 IY2 = IY 

465 INC = 6 
IX = IX+1 
GO TO 469 

466 INC = 7 
IY = IY+1 
IY2 = IY 
GO TO 458 

467 INC = 8 
IX = IX-1 

IF (IX) 468*468*469 
DID WE GO BACK TO STARTING POINT... 

468 IF ( ABS(XD-XS)-2. ) 481,481*490 

481 IF ( ABS( YD-YS)-2. ) 482*482*490 

C IF BACK AT START, PRINT STARTING COORD. 

482 WRITE (7,497) NSL,XS*YS 
GO TO 490 

470 GO TO (472,472,473,474,475*476,477,468), INC 

472 IY = IY+1 
INC = 3 

471 IY2 = IY 
NS = -1 
GO TO 420 

473 IY = IY-2 
INC = 4 
GO TO 471 

474 INC = 5 
1X2 = IX 

479 IY2 = IY+1 

480 NS = 1 

GO TO 420 

475 INC = 6 
IY = IY+1 
GO TO 479 

476 INC = 7 
IX = IX+1 
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non 



1X2 = IX 
GO TO 480 
477 INC = 8 

IY = IY-1 
GO TO 479 

END OF SHORELINE SEGMENT 
490 WRITE (7,499) 

499 FORMAT <5X, 22H END SHORELINE SEGMENT ) 
C GO LOOK FOR ANOTHER SEGMENT 
GO TO 404 
495 RETURN 
END 
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APPENDIX (C) 



FORMAT FOR COMMAND /DATA INPUT 



CARD COLUMNS 

1 20 21-25 26-30 31-40 41-50 



ANALYTICAL MATRIX 1X1 

I V 1 

LIMITS OF TOPO I XI 

I Y 1 

TOPO IN FEET FOLLOWS 
TOPO IN FATH FOLLOWS 

I NCREMFNT 

RAY DATA 1 NTC 

RAY DATA 2 

WATER LEVEL 

TRACE RAY 

LOCATE SHORELINE 

BREAKWATER 

BREAKWATER 

GAP 

SIZE GR 

WAVE LENGHTS DIFFR 
t.FNGTH DIFFRACTED 
ELIMINATF brkwtrs 



1X2 DX DM 

I Y 2 DY 

1X2 
I Y 2 

S 

T AD 

81 BN 

FI 

BX AS 8 YAS 

BXCS BYCS 



WLD 

WLD 



ID 



I XN 



IYN 



GRID 



DIFFRACT A 
DIFFRACT B 
DIFFRACT C 
DIFFRACT D 
DIFFRACT GAP 

DIFFR COORD X Y 



51-60 



X 



BX BS 
BXDS 



(CONTINUED) 



61-70 



Y 



BYBS 

BYDS 
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1 20 


CARD COLUMNS 
21-25 26-30 31-40 


41-50 


51-60 


61-70 


WAVE SPEED 


C 








CALC WAVE SPEED 


X 


Y 






( BLANK ) 


N 1 N2 F 1 


F 2 


F 3 


F 4 


all DONE 











NOTE. ..RIGHT JUSTIFY INTEGER DATA ITEMS 
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APPENDIX (D) 



REFRACTION TEST ON ANALYTICAL BEACH 
(SEE TABLF 1. ) 



ANALYTICAL DATA FROM X= 1 TO 50 AND Y= 1 TO 40 
T FST CASE 1 PERIOD = 4.0SEC. ANGL E= 45.0DEG.X= 3.00Y= 3.00 



N 


X 


Y 


A 


C 




BETA 


COEF 

RFTN 


shoal 

COEF 


CUR V 


S 


DPTH 


1 


8.000 


i.ooo 


45.000 


20.49 


1 


.000 


1 .000 


1 .00 


0.0001 


4.00 


56.00 


2 


5.828 


5.829 


45.019 


20.48 


1 


.000 


1 .000 


1.00 


0.0002 


4. CO 


50.34 


3 


7.241 


7.244 


45.040 


20.47 


1 


.00] 


1 .000 


1.00 


O.C003 


2.00 


47.51 


4 


8.654 


8.659 


45.074 


20.46 


1 


.001 


0.999 


1.00 


0.0004 


2.00 


44.68 


5 


10.066 


10.076 


45.125 


20.43 


1 


.002 


0.999 


0.99 


0.0007 


2.00 


41.85 


6 


1 1 .476 


1 1 .495 


45.222 


20.40 


1 


.00a 


0.998 


0.99 


0.0010 


2.00 


39.01 


7 


12.883 


12.916 


45.369 


20.35 


1 


.006 


0.997 


0.98 


0.0016 


2.00 


36.17 


8 


14.285 


14.342 


45.^8^ 


20.27 


1 


.010 


0.995 


0.98 


0.0023 


2.00 


33.32 


9 


15.681 


15.775 


45.917 


20.15 


1 


.016 


0.°92 


0.97 


0.0035 


2.00 


30.45 


10 


17.066 


17.217 


46.402 


19.97 


1 


.024 


0.988 


0.96 


0.0051 


2.00 


27.57 


1 1 


17.753 


17.943 


46.728 


19.85 


1 


.029 


0.986 


0.96 


0.0061 


1.00 


26.11 


12 


18.436 


18.674 


47.1 ] 5 


19.71 


1 


.036 


0.983 


0.95 


0.0074 


1.00 


24.65 


1 3 


19.114 


] 9. 409 


47.579 


19.54 


] 


.041 


0.979 


0 . 94 


0.0089 


1.00 


23.18 


14 


19.785 


20.151 


48 • 1 1 a. 


19.33 


1 


.052 


0.975 


0.94 


0.0106 


1.00 


21.70 


15 


20.118 


20.524 


48.451 


19.21 


1 


.058 


0.972 


0.94 


0.0115 


0.50 


20.95 


16 


20.448 


20.899 


48.798 


19.08 


1 


.063 


0.970 


0.93 


0.0126 


0.50 


20.20 


17 


20.776 


21.277 


49.175 


18.93 


1 


.069 


0.967 


0.93 


0.0137 


0.50 


19.45 


18 


21.102 


21.656 


49.586 


18.78 


1 


.076 


0.964 


0.93 


0.0150 


0.50 


18.69 


19 


21.424 


22.038 


50.033 


18.60 


1 


.083 


0.961 


0.92 


0.0163 


0.50 


17.92 


20 


21 .744 


22.423 


50.520 


18.41 


1 


.091 


0.957 


0.92 


0.0177 


0.50 


17.15 


2 1 


22.060 


22.810 


51.050 


18.21 


1 


.099 


0.954 


0.92 


0.0193 


0.50 


16.38 


22 


22.^72 


23.200 


51.627 


17.98 


1 


. 108 


0.^50 


0.92 


0.0210 


0.50 


15.60 


23 


22.68] 


23.594 


5 2 • 2 5 F 


17.73 


1 


. 1 17 


0.946 


0.92 


0.0229 


0.50 


14.81 


24 


22.984 


21 . 991 


52.939 


17.45 


1 


.128 


0.942 


0.91 


0.0249 


0.50 


14.02 


25 


23.283 


24.392 


51.683 


17.15 


] 


.139 


0.937 


0.91 


0.0271 


0.50 


13.22 


26 


23.576 


24.797 


54.495 


16.82 


] 


.150 


0.912 


0.91 


0.0296 


0.50 


12.41 


27 


23.864 


25.206 


55.379 


16.45 


1 


.163 


0.927 


0.92 


0.0323 


0.50 


11.59 


28 


24 . 144 


25.620 


56.346 


16.05 


] 


.176 


0.922 


0.92 


0.0353 


0.50 


10.76 


29 


24.417 


26.039 


57.401 


15.60 


1 


.191 


0.916 


0.92 


0.0386 


0.50 


9.92 


30 


24.682 


26.463 


58.563 


15.11 


1 


.206 


0.911 


0.93 


0.0424 


0.50 


9.07 


3 1 


24.938 


26.892 


59.84] 


14.55 


1 


.222 


0.905 


0.93 


0.0468 


0.50 


8.22 


32 


25. 184 


27.328 


61.253 


13.93 


1 


.239 


0.898 


0.94 


0.0519 


0.50 


7.34 


33 


25.419 


27.769 


62.824 


13.23 


1 


.257 


0.892 


0.96 


0.0579 


0.50 


6.46 


34 


25.640 


28.218 


64.586 


12.43 


1 


.276 


0.885 


0.98 


0.0653 


0.50 


5.56 


35 


25.847 


28.671 


66.585 


11 .51 


1 


.297 


0.878 


1.00 


0.0745 


0.50 


4.65 


36 


26.036 


29. 136 


68.886 


10.43 


1 


.118 


0.871 


1.04 


0.0869 


0.50 


3.73 


37 


26.205 


29.606 


71.610 


9.14 


1 


.341 


0.864 


1.10 


0.1048 


0.50 


2.79 


38 


26.349 


30.085 


74.994 


7.50 


1 


.365 


0.856 


1.20 


0.1347 


0.50 


1.83 


39 


26.458 


30.573 


79.70Q 


5.19 


] 


. 390 


0.848 


1.42 


0.2045 


0.50 


0.85 



RAY STOPPED 



-lllS- 



APPENDIX (E) 



SAMPLE PROGRAM RUN 



DIFFRACTION OPERATIONS CALCULATE VALUES FOR AN ASSORTMENT OF 
MISCELLANEOUS CASES. 



REFRACTION OPERATIONS 


SHOW 


VIRGIN I A 


BEACH TESTS 


(SEE FIG. 


4 ) 


T NPlJT DATA OFf K 












WAVE SPEED 




20. 








SIZF GRID 


50 


40 80. 








BW 




10. 


10. 


20. 


20. 


BW 




22. 


22. 


30. 


30. 


RAY DATA 1 
DIFFRACT A 


9 


4 • 


90. 


1 . 


1 . 


DIFFR COORD 




8. 


14. 










10. 


13. 










12. 


14 . 






RAY DATA 1 


11 


4 . 


135. 


1 . 


1. 


DIFFR COORD 




6 . 


12. 










6.5 


12.5 










7. 


13. 










8. 


14. 










12. 


13. 






DIFFRACT D 
DIFFR COORD 




25. 


30 . 










26. 


32. 










27. 


33. 










27.5 


33.5 










28. 


34. 






RAY DATA 1 
DIFFRACT GAP 


1 


4.0 


90. 


1 . 


1. 


DIFFR COORD 




21. 


26. 










23. 


25 . 










20.5 


23. 






ELIMIN BREAKWATERS 
BW 




10. 


10 . 


20. 


10. 


DIFFRACT A 
RAY DATA 1 


1 


4. 


90. 


2. 


2. 


DIFFR COORD 




9.9 


15 . 










10. 


15 . 










10.1 


15. 






DIFFRACT B 
DIFFR COORD 




19.9 


15 . 










20. 


15 . 










20.1 


15. 







eliminate BRKWTRS 
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SIZE GRID 50 

ANALYTICAL MATRIX 1 

33 

LIMITS OF TOPO 1 

1 

TOPO IN FATH FOLLOWS 



9 5 


95 


96 


97 


97 


92 


90 


88 


85 


83 


82 


82 


82 


82 


82 


83 


83 


84 


85 


85 


88 


89 


98 


92 


92 


92 


93 


94 


95 


95 


92 


89 


87 


84 


81 


81 


81 


81 


80 


80 



40 3000. 

50 0. -24. 

40 0 . 

50 

32 



98 


98 


99 


96 


95 


82 


82 


82 


82 


8? 


81 


82 


82 


82 


82 


86 


86 


87 


87 


88 


93 


93 


94 


93 


92 


96 


96 


97 


96 


94 


80 


80 


80 


81 


81 


80 


80 


80 


81 


81 



( INTFRIM CARDS OMITTED) 



-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-4 0 


-40 


-40 


-40 


-40 


-40 


-40 


435 


44 


44 5 


45 


452 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-35 


-30 


-25 


12 


19 


22 


28 


30 


405 


41 


415 


42 


425 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-1 5 


-10 


-00 


08 


12 



-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-35 


-30 


-25 


-20 


-15 


455 


458 


46 


465 


47 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-20 


-15 


-10 


-00 


09 


33 


36 


38 


40 


405 


435 


435 


44 


445 


45 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-40 


-35 


-30 


-25 


-20 


19 


20 


22 


25 


29 



(INTERIM CARDS OMITTED) 



-40 -40 


-40 


-40 


-40 


-40 


-40 -40 


-40 


-40 


-40 


-40 


. -40 -40 


-40 


-40 


-40 


-40 


-40 -40 


-40 


-40 


-40 


-40 


-40 -40 


-40 


-40 


-40 


-40 


RAY DATA 1 






31 


4* 


TRACE RAY 
RAY DATA 1 






32 


4. 


TRACE RAY 
RAY DATA 1 
TRACE RAY 






33 


4. 



LOCATE SHORELINE 
ALL DONE 



-40 


-40 


-40 


-40 




-40 


-40 


-40 


-40 




-40 


-40 


-40 


-40 




-40 


-40 


-40 


-40 




-40 


-40 


-30 


-28 






30. 


5.0 




11.1 




30 . 


6.0 




9.3 




30. 


7.0 




7.9 



-If*)- 



. 






V 



OUTPUT INFORMATION DECK 



WAVE SP^ED IS SFT AT 20.00 FT/S 

GRID IS 50 PY 40 WITH 80.00 FOOT SQUARES. 
RRFAKWAT f R FROM 1°.00, 10.00 TO 20.00. 20.00 
BREAKWATER FROM 22.00, 22.00 TO BO. 00, 30.00 
DIFFRACTION CALCULATIONS FOR A END OF BREAKWATER. 
WAVE CHARACTERISTICS 



SPEED20.00 FT/SEC. 



ANGLE 90.00 DFG • , PFRIOD 4.00 SEC., 

AT POINT X = 8.00 Y = 14.00 

DIR OF TRAVEL IS 92.37 DEG. ,COEF OF DIFFRACTION IS 1.1040 
WAVE CHARACTERISTICS 

ANGLE 90.00 DFG., PFRIOD 4.00 SFC.» SPEED20.00 FT/SEC. 

AT POINT X =10.00 Y= 13.00 

DIR OF TRAVFL IS 82.60 DEG. , COEF OF DIFFRACTION IS 0.5490 
WAVE CHARACTERISTICS 

ANGLE 90.00 DEG., PERIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X =12.00 Y= 14.00 

DIR OF TRAVEL IS 62.67 DEG., COEF OF DIFFRACTION IS 0.2278 
WAVE CHARACTERISTICS 

ANGLE 135.00 DEG., PERIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X = 6.00 Y= 12.00 

DIR OF TRAVEL IS 137.41 DEG., COEF OF DIFFRACTION IS 1.0758 
WAVE CHARACTERISTICS 

ANGLE 135.00 DFG • > PFRIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X = 6.50 Y= 12.50 

DIR OF TRAVEL IS 134.01 DEG., COEF OF DIFFRACTION IS 0.8031 
WAVE CHARACTERISTICS 

ANGLE 135.00 DEG., PERIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X = 7.00 Y= 13.00 

DIR OF TRAVEL IS 128.80 DEG., COEF OF DIFFRACTION IS 0.5283 
WAVE CHARACTERISTICS 

ANGLE 135.00 DEG., PFRIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X = 8.00 Y= 14.00 

DIR OF TRAVEL IS 115.00 DFG.,COEF OF DIFFRACTION IS 0.2466 
WAVE CHARACTERISTICS 

ANGLF 135.00 DFG., PERIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X =12.00 Y= 13.00 

DIR OF TRAVEL IS 56.31 DEG., COEF OF DIFFRACTION IS 0.1203 
DIFFRACTION CALCULATIONS FOR D END OF BREAKWATER. 

WAVE CHARACTERISTICS 

ANGLE 135.00 D C G . » PERIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X =25.00 Y= 30.00 

DIR OF TRAVEL IS 180.17 DEG., COEF OF DIFFRACTION IS 0.1309 
WAVE CHARACTERISTICS 

ANGLF 135.00 DFG., PERIOD 4.00 SEC., SPEED20.00 ET/SEC. 

AT POINT X =26.00 Y= 32.00 

DIR OF TRAVEL IS 155.00 DFG., COEF OF DIFFRACTION IS 0.2466 
W A v c characteristics 
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ANGLF 135.00 DFG • , PFRIOO 4.00 SFC.» SPEFD20.00 FT/SEC. 

AT POINT X =27.00 Y= 33.00 

DIR OF TRAVEL IS 141.20 DFG., COEF OF DIFFRACTION IS 0.5263 
WAVE CHARACTERISTICS 

A NGL r 135.00 D P G . , PERIOD 4.00 SEC.. SPEED20.00 FT/SEC. 

AT POINT X =27.50 Y= 33.50 

DIR OF TRAVEL IS 135.99 DEG..COEF OF DIFFRACTION IS 0.8031 
WAVE characteristics 

ANGLF 135.00 DFG.» PERIOD 4.00 SEC.. SPEED20.00 FT/SEC. 

AT POINT X =28.00 Y= 34.05 

DIR OF TRAVEL- IS 132.59 DEG..COEF OF DIFFRACTION IS 1.0758 
DIFFRACTION CALCULATIONS FOR BREAKWATER GAP. 

BASED ON IMAGINARY GAP ASSUMED ^ORmAL TO WAVE RAY. 

IMAGINARY GAP FROM 20.000,21.000 TO 22.000,21.000 

WAVE characteristics 

ANGLF 90.00 DFG • » PFRIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X =21.00 Y= 26.00 

DIR OF TRAVEL IS 90.00 DFG.,COEF OF DIFFRACTION IS 0.8493 
WAVE CHARACTERISTICS 

ANGLF 90.00 DFG., PERIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X =23.00 Y= 25.00 

DIR OF TRAVEL IS 77.82 D C G • , COFF OF DIFFRACTION IS 0.2386 
WAVE CHARACTERISTICS 

ANGLF 90.00 DFG., PFRIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X =20.50 Y= 23.00 

DIR OF TRAVEL IS 94.25 DEG. , COEF OF DIFFRACTION IS 0.8883 
ELIMINATED BREAKWATERS 

BREAKWATER FROM 10.00, 10.00 TO 20.00, 10.00 

DIFFRACTION CALCULATIONS FOR A END OF BREAKWATER. 

WAVE CHARACTERISTICS 

ANGLE 90.00 DEG., PERIOD 4.00 SEC., SPEED20.00 FT/SEC. 

AT POINT X = 9.90 Y= 15.00 

DIR OF TRAVEL IS 84.98 DEG., COEF OF DIFFRACTION IS 0.5569 
WAVE CHARACTERISTICS 

ANGLF 90.00 DFG., PFRIOD 4.00 SEC., SPFFD20.00 FT/SFC. 

AT POINT X =10.00 Y = 15.00 

DIR OF TRAVEL IS 84.28 DEG., COEF OF DIFFRACTION IS 0.5260 
WAVE CHARACTERISTICS 

ANGLE 90.00 DFG. , PERIOD 4.30 SEC., SPEED20.00 FT/SEC. 

AT POINT X =10.10 Y = 15.00 

DIR OF TRAVEL IS 83.56 DEG., COEF OF DIFFRACTION IS 0.4968 
DIFFRACTION CALCULATIONS FOR B END OF BREAKWATER. 

WAVE CHARACTERISTICS 

ANGLF 90.00 DFG., PERIOD 4.00 SFC., SPEED20.00 FT/SFC. 

AT POINT X =19.50 Y = 15.00 

DIR OF TRAVEL IS 96.44 DFG.,COFF OF DIFFRACTION IS 0.4968 
WAVF CHARACTERISTICS 

ANGLF 90.00 DFG., PFRIOD 4.00 S C C., SPEED20.00 FT/SEC. 

AT POINT X =20.50 Y= 15.50 

DIR OF TRAVEL IS 95.72 DFG.,COEF OF DIFFRACTION IS 0.5260 
WAVE CHARACTERISTICS 
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ANGLE 90.00 DFG.» PERIOD A. 00 SEC.* SPEED20.00 FT/SEC. 

AT POINT X =20.10 Y= 15.00 

DIR OF TRAVEL IS 95. °2 DFG. * COEF OF DIFFRACTION IS 0.5569 
FLr-'INATED prfakwaters 

GRID IS 50 PY 40 WITH 300C.00 FOOT SQUARES. 

ANALYTICAL DATA FROM X= 1 TO 50 AND Y= 33 TO 40 
TOPO DATA FROM X= 1 TO 50 AND Y= 1 TO 32 
TEST CASE 31 PERIOD = 4.0SEC. ANGL E= 30.0DEG.X= 5.00Y=11.10 



N 


X 


Y 


A 


C 


BETA 


COEF 

RETN 


shoal 

COEF 


CUR V 


S 


DPTH 


1 


5 .000 


11.100 


30.000 


20.42 


1 .COO 


1.000 


1.00 


0.0014 


4.00 


40.64 


2 


6.731 


12.102 


30.159 


2 0.43 


1 .009 


0.995 


1 . 00 


0.0012 


2.00 


41.82 


-x 


8.459 


13. i 09 


30.2*5 1 


2 0.44 


1.024 


0.Q88 


1.00 


0.0009 


2.00 


42.10 


4 


10.185 


14.119 


30.392 


2 0.43 


1 .045 


0.Q78 


1.00 


C.0017 


2.00 


41.25 


5 


1 1 .908 


15.135 


30.622 


20.39 


1 .0 76 


0.964 


1.00 


0.0023 


2.00 


3 8.54 


6 


13 .626 


16.158 


30.923 


20.36 


1.112 


0.948 


0.99 


0.0029 


2.00 


36.61 


7 


15.539 


17.191 


31.28? 


20. ^0 


1.155 


0 . ° 3 1 


0.99 


0.0033 


2.00 


34.56 


8 


17.044 


18.236 


31 .683 


20.26 


1 .202 


0.912 


0.99 


0.0036 


2.00 


33.17 


9 


18.742 


19.2*52 


32.110 


20.21 


1.261 


0.891 


0.98 


0.0039 


2.00 


31. 8 Q 


10 


20.429 


20.367 


32.875 


20.09 


1 .347 


0.862 


0.98 


0.0113 


2.00 


29.43 


1 1 


20.848 


20.640 


33.244 


20.06 


1.382 


0.851 


0.97 


0.012° 


0.50 


28.94 


12 


21.266 


20.01 5 


33.653 


20.04 


1 .425 


0.838 


0.97 


0.0160 


0.50 


28.57 


13 


21.681 


21 . 194 


34.177 


19.99 


1.476 


0.823 


0.97 


0.0211 


0.5 0 


27.87 


14 


22.092 


21.477 


34.871 


19. *50 


1 .53*5 


0.806 


0.97 


0.1*277 


0. 50 


2 6.64 


1 5 


22.500 


21.767 


35.786 


19.74 


1.618 


0.786 


0.96 


0. J368 


0.5C 


24.98 


16 


22.903 


22.063 


37.001 


19.55 


1.716 


0.763 


0. 05 


0.0487 


0.50 


2 3.31 


1 7 


23.298 


22.370 


38.622 


1 9.^4 


1 .840 


0.737 


0.95 


0.0657 


0.50 


21.79 


1 8 


23.683 


22.68*5 


40.761 


19.0*5 


2.00? 


0.707 


0.94 


0.0840 


0.50 


20.30 


19 


24.053 


23.0?= 


43.608 


18.75 


2.216 


0.672 


0.93 


0. 1087 


0.50 


18.57 


20 


24.405 


23.380 


47.1 37 


18.25 


2.504 


0.632 


0.93 


0. 1493 


0.50 


16.55 


21 


24.728 


23.762 


52.329 


1 7.37 


2.905 


0.587 


0.02 


0.2210 


0.50 


13.73 


22 


25.004 


24. 179 


60.757 


15.41 


3.48° 


0.535 


0.93 


0.3924 


0. 50 


9.59 


23 


25.169 


24.650 


80.545 


9.50 


4 . 397 


0.477 


1.09 


1.1307 


0.50 


3.10 


RAY STOPPED 
TEST CASE 32 P c 


■RIOD = 


4.0SEC. ANGL E= 30.0DEG 


. X = 6 • OOY= 9. 


,30 


N 


X 


Y 


A 


C 


BETA 


COEF 

RFTN 


shoal 

COEF 


CUR V 


S 


DPTH 


1 


6.000 


9. 300 


30.000 


20.45 


1.000 


1 . OQ 0 


1.00 


0.0004 


4.0 0 


44.06 


2 


7.7^2 


10.301 


30.051 


2 0.45 


1 .000 


1.000 


1.00 


0.0005 


2.00 


4 3*94 


8 


*5.462 


1 1 . 303 


30.104 


2 0.44 


1 .000 


1.000 


1.00 


0. 0002 


2.00 


42.85 


4 


11.192 


12.307 


’0.12’ 


20.44 


1.000 


1 .000 


1.00 


0.0001 


2.00 


42.52 


5 


12.922 


13.3]0 


30.1 37 


20.44 


1.000 


1 .00*5 


1.00 


0.0000 


2.00 


42.64 


6 


14.652 


14.315 


30.1 40 


20.44 


1 .00? 


0 • °99 


1.00 


O.C007 


2.00 


42.15 


7 


16.380 


15.321 


30.264 


2 0.41 


1.008 


0 . ° 96 


1.00 


0.0014 


2.00 


40.02 


8 


18.106 


16.332 


30.471 


20.37 


1.018 


0.991 


0.99 


0.0022 


2.00 


3 7.15 


9 


19.828 


17.350 


30. 706 


20.31 


1.030 


0.985 


0.99 


0.0015 


2.00 


34.91 


10 


21.546 


1 R • ^ 7 3 


30.867 


20.28 


1.040 


0.081 


0.98 


0.0015 


2.00 


3 Q • 79 


1 1 


23.260 


19.403 


31.105 


20 . 2 ? 


1.052 


0.°75 


0.98 


0.0029 


2.00 


32.3 ] 


1 2 


24.970 


20.441 


31.4^2 


2 0.16 


1.06Q 


0.967 


0.98 


0.0022 


2.00 


31.05 


1 3 


26.673 


21.490 


31 .844 


20.08 


1.091 


0.957 


0.97 


0.0061 


2.00 


29.28 


1 4 


27.520 


22.021 


32.245 


20.04 


1 . I 1 1 


O .049 


0.97 


0.0071 


1 .00 


28.51 
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1 5 


28.364 


22.557 


32.667 


19.Q9 




. 1 3 7 


0.938 


0.97 


0.0074 


1 . 00 


2 7.91 


1 6 


29.203 


23.101 


33.701 


19.01 


1 


.176 


O.Q?? 


0.96 


0.0126 


1.00 


26.76 


17 


29.621 


23.376 


33.^86 


19.86 


1 


. 203 


0.912 


0.96 


0.0135 


0.50 


26.2 5 


1 8 


30. 036 


23.654 


33.090 


19.86 


1 


.236 


0.899 


0.96 


0.0147 


0.50 


25.98 


19 


30.45C 


23.936 


34.478 


19.77 


1 


.277 


0.885 


0.96 


0.0204 


0.50 


2^.29 


20 


30.860 


24.221 


35.200 


19.63 


1 


.327 


0.868 


0.95 


0.0313 


0. 50 


23.92 


21 


31 . 266 


24.513 


36.293 


19.41 


1 


. 391 


0.848 


0.95 


0.0458 


0.50 


22.23 


22 


31 .665 


24.815 


37.846 


19.08 


1 


.474 


0.826 


0.94 


0.0633 


0.50 


2 0.21 


23 


32.054 


25.128 


39.858 


18.72 


1 


.581 


0.795 


0.93 


0.0758 


0. 50 


18.42 


24 


32.431 


25.687 


42.276 


18.36 


1 


.7 18 


0.763 


0.93 


0 .0944 


0. 50 


16.96 


25 


32.791 


25.804 


46,664 


17.6° 


1 


.896 


0.726 


0.92 


0.1501 


0.50 


It .70 


26 


33.122 


26. 179 


51,534 


16.45 


2 


. 14° 


0.682 


0.92 


0.2796 


0.5C 


11.59 


27 


3 •a . 393 


26.599 


62.640 


14.04 


2 


.574 


0.623 


0.95 


0.5250 


0. 50 


7.49 


28 


33.514 


27.084 


89.359 


6.89 


3 


.34? 


0.567 


1.25 


1.5308 


0. 50 


1.5? 


RAY 


STOPPED 




















TEST CASE 


33 PERIOD = 


4.0SEC 


. ANGLE* 30.0DEG 


• X= 7 . 00Y= 7. 


,90 
















COE F 


SHOAL 








N 


X 


Y 


A 


C 




BETA 


R FT N 


COE F 


CUR V 


S 


DPTH 


1 


7.300 


7.900 


39.000 


20.47 


1 


.000 


1.000 


1.00 


0.0003 


4.00 


46.95 


2 


8.732 


8.901 


30.037 


20.67 


1 


.001 


1.000 


1.00 


0.0003 


2.00 


46.19 


3 


10.463 


9. Q02 


30.066 


20.46 


1 


.00? 


0.999 


1.00 


0.0C04 


2.00 


44.79 


4 


12.193 


10.905 


30.108 


2 0.44 


1 


.004 


0.998 


1.00 


0.0003 


2.00 


42.94 


5 


1^.923 


11 .908 


30.137 


20.64 


1 


.005 


0.997 


1.00 


-0.0000 


2.00 


42.04 


6 


15.653 


12.913 


30.136 


2 0.63 


1 


. 1 06 


0.Q97 


1.00 


-9.0002 


2.00 


41.96 


7 


17.383 


13.016 


30.115 


2 0.44 


1 


.006 


O .997 


1.00 


0.0001 


2.00 


42.63 


8 


19.112 


14.Q21 


30.176 


20.42 


1 


.0 09 


0.995 


1.00 


0.0011 


2.0C 


40.99 


9 


20.840 


15.929 


30.359 


20.38 


1 


. H6 


0.992 


0.99 


0.0021 


2.00 


38.00 


10 


22.563 


16.943 


30.504 


20.37 


1 


.024 


0.988 


0.99 


0.0019 


2.00 


35.3? 


1 1 


24.283 


17.964 


30.790 


20.30 


1 


.029 


0.986 


0.98 


0.0012 


2. 00 


94.23 


12 


26.000 


18.991 


30.954 


20.26 


1 


.038 


0.983 


0.98 


0.0020 


2.00 


33.01 


1 3 


27.712 


20.023 


31.214 


20.21 


1 


.046 


0.978 


0.98 


0.0024 


2.00 


31.78 


14 


29.42C 


21.064 


31.525 


20.17 


1 


.060 


0.971 


0.97 


0.0031 


2.00 


30.21 


1 5 


31.121 


22.116 


31.912 


20.06 


1 


.077 


0 . Q 6 3 


0.97 


0.0035 


2.00 


28.95 


16 


32.815 


23.180 


32 .389 


19.97 


1 


.101 


O .933 


0.96 


0.0051 


2.00 


27.60 


1 7 


33.657 


23.718 


32.740 


19.92 


1 


.118 


0.946 


0.96 


0.0079 


1.00 


26.93 


18 


34.496 


24.262 


33.175 


19.86 


1 


. 143 


0.935 


0.96 


0.007? 


1 . 00 


26.16 


19 


35.33] 


24.814 


33.726 


19.78 


1 


. 1 78 


0.922 


0.96 


0.0135 


1.00 


25.37 


20 


35.746 


25.CQ3 


34.205 


1°.69 


1 


.203 


0.912 


0.95 


0.0200 


0. 50 


24.47 


21 


36.157 


25.376 


34.868 


19.58 


1 


.236 


0.899 


0.95 


0.0263 


0. 50 


23.49 


22 


36.565 


25.665 


35.720 


19.47 


1 


.280 


0.884 


0.95 


0.0334 


0.50 


22.72 


23 


36.969 


25.961 


36.864 


19.34 


1 


.339 


0.864 


0.94 


0.0465 


0.50 


21.80 


24 


37.364 


26.267 


38.484 


19.06 


1 


.41° 


0.839 


0.93 


0.0706 


0.5 0 


20.0° 


25 


37.749 


26.586 


41.052 


18.50 


1 


. 53? 


0.808 


0.93 


0.1130 


0.50 


17.51 


26 


00 

• 

h-J 


26.928 


45.200 


17.52 


1 


.698 


0.768 


0.92 


0. 1835 


0.50 


14.20 


27 


38.444 


27.304 


52.1 73 


15.81 


1 


. 0 44 


0.717 


0.92 


0.3198 


0.50 


10.30 


28 


38.703 


27.731 


65. ^06 


12.4 1 


2 


. 346 


0.653 


0.98 


0.6438 


0.5 0 


5. 5 C 


RAY 


STOPPED 






















SHORELINE COORDINATES 
















1 








1.000 


16.048 












2 








2.000 


16. 500 












3 








3.000 


17.000 
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4 


3.00C 


17.000 


S 


4.000 


17 .464 


6 


5.000 


17.667 


7 


5.588 


18.000 


8 


6.000 


18.318 


9 


7.000 


18 . 600 


10 


8.000 


19.000 


11 


8.000 


19.000 


1? 


9.000 


1 9.444 


13 


10.000 


19.643 


14 


11.000 


20.000 


15 


11.000 


20.000 


16 


12.000 


20.400 


17 


13.000 


20.783 


18 


14.000 


21.000 


19 


14.000 


21 . 000 




(INTERIM CARDS OMITTED) 




52 


40.000 


28.400 


53 


41 .000 


28.520 


54 


42 . 000 


28.643 


55 


43.000 


28.800 


56 


44.000 


29.000 


57 


44.000 


29.000 


58 


45.000 


29.400 


69 


46.000 


29.556 


60 


47.000 


29.667 


61 


48.000 


29.815 


62 


49.000 


30.000 


66 


49.000 


30.000 


64 


50.000 


30.400 



END SHORFLINF SEGMENT 



ALL DONF. 
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